1 Introducción y definición de objetivos

El análisis de datos en el mercado inmobiliario es algo muy común y que lleva muchísimos años desarrollándose. Con el objetivo de predecir que aspectos influyen principalmente en el precio de las casas en Madrid, hemos seleccionado un dataset (kaggle - Madrid House Price) con viviendas en venta de la capital española.

library(readr)
library(dplyr)
library(tidyr)
library(kableExtra)
library(ggplot2)
library(GGally)
library(gridExtra)
library(cowplot)
library(ggcorrplot)
library(gmodels)
library(caret)
library(ggfortify)
library(scales)
library(class)
library(distances)
library(visreg)
library(rpart)
library(rpart.plot)
library(rattle)
library(randomForest)
library(e1071)
library(pROC)
library(cluster)
library(tidyverse)
library(stringr)
library(Metrics)
library(factoextra)
library(NbClust)
library(neuralnet)

# Librerías para SERIES TEMPORALES
library(xts)
library(forecast)
library(tsibble)
library(tidyquant)
library(feasts)
library(fable)
library(stats)
library(zoo)
library(DT)
paste(R.Version()$version.string)
## [1] "R version 4.1.2 (2021-11-01)"

2 Análisis exploratorio inicial

2.1 Lectura y preparación de los datos

En esta primera etapa cargamos el dataset descargado y hacemos una observación superficial de las variables recogidas y sus características. Se representan también las 10 filas iniciales.

mhp <- read_csv("house_price_madrid_14_08_2022.csv")
head(mhp, 10) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
price house_type house_type_2 rooms m2 elevator garage neighborhood district
495000 planta 1 exterior 3 118 TRUE TRUE Chopera Arganzuela
485000 planta 2 exterior 2 82 TRUE TRUE Palos de Moguer Arganzuela
315000 planta 2 exterior 2 72 FALSE FALSE Legazpi Arganzuela
585000 planta 4 exterior 2 174 TRUE TRUE Palos de Moguer Arganzuela
255000 bajo exterior 3 75 FALSE FALSE Acacias Arganzuela
299000 planta 1 exterior 1 69 TRUE FALSE Chopera Arganzuela
265000 planta 3 exterior 2 54 TRUE FALSE Delicias Arganzuela
290000 planta 3 interior 4 69 TRUE FALSE Chopera Arganzuela
660220 planta 2 exterior 3 129 TRUE TRUE Imperial Arganzuela
525000 planta 4 exterior 4 111 TRUE FALSE Chopera Arganzuela
summary(mhp)
##      price           house_type        house_type_2           rooms       
##  Min.   :     725   Length:15975       Length:15975       Min.   : 1.000  
##  1st Qu.:  195000   Class :character   Class :character   1st Qu.: 2.000  
##  Median :  359973   Mode  :character   Mode  :character   Median : 3.000  
##  Mean   :  624233                                         Mean   : 2.848  
##  3rd Qu.:  749000                                         3rd Qu.: 3.000  
##  Max.   :13950000                                         Max.   :41.000  
##        m2         elevator         garage        neighborhood      
##  Min.   :  1.0   Mode :logical   Mode :logical   Length:15975      
##  1st Qu.: 66.0   FALSE:4597      FALSE:11528     Class :character  
##  Median : 93.0   TRUE :11378     TRUE :4447      Mode  :character  
##  Mean   :124.8                                                     
##  3rd Qu.:142.0                                                     
##  Max.   :989.0                                                     
##    district        
##  Length:15975      
##  Class :character  
##  Mode  :character  
##                    
##                    
## 
str(mhp) 
## spc_tbl_ [15,975 × 9] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ price       : num [1:15975] 495000 485000 315000 585000 255000 ...
##  $ house_type  : chr [1:15975] "planta 1" "planta 2" "planta 2" "planta 4" ...
##  $ house_type_2: chr [1:15975] "exterior" "exterior" "exterior" "exterior" ...
##  $ rooms       : num [1:15975] 3 2 2 2 3 1 2 4 3 4 ...
##  $ m2          : num [1:15975] 118 82 72 174 75 69 54 69 129 111 ...
##  $ elevator    : logi [1:15975] TRUE TRUE FALSE TRUE FALSE TRUE ...
##  $ garage      : logi [1:15975] TRUE TRUE FALSE TRUE FALSE FALSE ...
##  $ neighborhood: chr [1:15975] "Chopera" "Palos de Moguer" "Legazpi" "Palos de Moguer" ...
##  $ district    : chr [1:15975] "Arganzuela" "Arganzuela" "Arganzuela" "Arganzuela" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   price = col_double(),
##   ..   house_type = col_character(),
##   ..   house_type_2 = col_character(),
##   ..   rooms = col_double(),
##   ..   m2 = col_double(),
##   ..   elevator = col_logical(),
##   ..   garage = col_logical(),
##   ..   neighborhood = col_character(),
##   ..   district = col_character()
##   .. )
##  - attr(*, "problems")=<externalptr>

El dataset contiene 15.975 observaciones (correspondientes a una vivienda cada una) y 9 variables (de las cuales son 6 cualitativas y 3 cuantitativas).

A continuación, la descripción de cada una de las variables:

  • price: precio

  • house_type: tipo de vivienda (casa, chalet, piso…)

  • house_type_2: si es exterior o interior

  • rooms: número de habitaciones

  • m2: metros cuadrados

  • elevator: si tiene ascensor

  • garage: si incluye garaje

  • neighborhood: barrio de Madrid

  • district: distrito de Madrid

Tras esta primera observación se harán algunas modificaciones en el dataset para transformar a variables de tipo categórico las que correspondan, así como transformar a binarias (1/0) aquellas con sólo dos valores posibles.

# Preparación de los datos. Convertimos a factores las variables que son categóricas y a 1-0 las binarias.

mhp <- mhp %>% rename(exterior = house_type_2)
mhp$exterior = ifelse(mhp$exterior == "exterior", 1, 0)
mhp$exterior = factor(mhp$exterior, levels = c(1,0))

mhp$elevator = ifelse(mhp$elevator == "TRUE", 1, 0)
mhp$elevator = factor(mhp$elevator, levels = c(1,0))

mhp$garage = ifelse(mhp$garage == "TRUE", 1, 0)
mhp$garage = factor(mhp$garage, levels = c(1,0))

mhp$house_type <- as.factor(mhp$house_type)
mhp$exterior <- as.factor(mhp$exterior)
mhp$elevator <- as.factor(mhp$elevator)
mhp$garage <- as.factor(mhp$garage)
mhp$neighborhood <- as.factor(mhp$neighborhood)
mhp$district <- as.factor(mhp$district)

2.2 Tratamiento de datos faltantes

Convertimos en NA aquellas observaciones que claramente contienen algún error:

  • La casa de 41 habitaciones (probablemente fueran 4 y esté mal imputado)

  • Las que tienen menos de 3 m2 (probablemente mal imputados al usar el punto para separar los miles)

  • La casa con un precio de 725 (claramente equivocado).

mhp$rooms[mhp$rooms > 40] <- NA
mhp$m2[mhp$m2 < 3] <- NA
mhp$price[mhp$price < 1000] <- NA

Por otro lado, la variable house_type_2 tiene 469 filas sin datos. Junto a los NA añadidos en el paso anterior son el 3.2% del total, por lo que optaremos por eliminar estas filas del dataset.

sum(is.na(mhp))/nrow(mhp)
## [1] 0.03205008
mhp <- na.omit(mhp)

2.3 División del dataset

Para finalizar la preparación de los datos, dividiremos nuestro dataset en tres partes, una de entrenamiento (train), con el que trabajaremos durante todo el análisis, otra de prueba (test) que nos servirá para comprobar la eficacia de los diferentes modelos predictivos que se apliquen, y finalmente uno de validación (validation), que se mantendrá sin usar hasta el último día del proyecto y con el cual validaremos el modelo final.

set.seed(108)
numero_total = nrow(mhp)
# Porcentajes de train, test y validation
w_train = .5
w_test = .25
w_validation = 1 - (w_train + w_test)

# Todos los índices
indices = seq(1:numero_total) 

# Muestreo
indices_train = sample(1:numero_total, numero_total * w_train)
indices_test = sample(indices[-indices_train], numero_total * w_test)
indices_validation = indices[-c(indices_train,indices_test)]

# Agrupamos

mhp_train = mhp[indices_train,]
mhp_test = mhp[indices_test,]
mhp_validation = mhp[indices_validation,]

3 Análisis univariante

3.1 Análisis de variables cualitativas

Una vez preparado el dataset comenzamos con el análisis de las variables cualitativas, que son las siguientes:

Tipo de vivienda

merge(setNames(as.data.frame(table(mhp_train$house_type)), c("house_type", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$house_type))*100, 2)), c("house_type", "prop (%)"))
) %>%
  arrange(desc(count)) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
house_type count prop (%)
planta 1 1581 20.45
planta 2 1332 17.23
planta 3 1180 15.26
bajo 996 12.88
planta 4 801 10.36
planta 5 545 7.05
planta 6 317 4.10
planta 7 196 2.53
chalet 180 2.33
entreplanta 131 1.69
planta 8 119 1.54
planta 9 78 1.01
casa 65 0.84
semi-sotano 48 0.62
planta 10 44 0.57
planta 12 23 0.30
planta 13 22 0.28
planta 11 21 0.27
sotano 19 0.25
planta 14 10 0.13
planta -1 8 0.10
planta 15 6 0.08
planta 16 4 0.05
planta 18 2 0.03
planta 19 2 0.03
planta 20 2 0.03
ggplot(mhp_train, aes(house_type)) +
  geom_bar(fill = "#0BB363") +
  coord_flip() +
  labs(x = "Tipo de vivienda", y = "Número de viviendas", title = "Viviendas por tipo") +
  theme(plot.title = element_text(hjust = 0.5))

Como era de esperar, las viviendas más habituales son pisos en plantas no demasiado altas (bajos, primeros, segundos…). El número de casas y chalets también es pequeño en proporción.

Orientación exterior

merge(setNames(as.data.frame(table(mhp_train$exterior)), c("exterior", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$exterior))*100, 2)), c("exterior", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
exterior count prop (%)
0 848 10.97
1 6884 89.03
ggplot(mhp_train, aes(exterior)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Exterior", y = "Número de vivienda de viviendas", title = "Viviendas exteriores") +
  theme(plot.title = element_text(hjust = 0.5))

Son una inmensa mayoría los pisos que dan a exterior.

Ascensor

merge(setNames(as.data.frame(table(mhp_train$elevator)), c("elevator", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$elevator))*100, 2)), c("elevator", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
elevator count prop (%)
0 2169 28.05
1 5563 71.95
ggplot(mhp_train, aes(elevator)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Ascensor", y = "Número de viviendas", title = "Viviendas con ascensor") +
  theme(plot.title = element_text(hjust = 0.5))

También es más frecuente que tengan ascensor, aunque casi el 30% no lo tienen.

Garaje

merge(setNames(as.data.frame(table(mhp_train$garage)), c("garage", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$garage))*100, 2)), c("garage", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
garage count prop (%)
0 5570 72.04
1 2162 27.96
ggplot(mhp_train, aes(garage)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Garaje", y = "Número de viviendas", title = "Viviendas con garaje") +
  theme(plot.title = element_text(hjust = 0.5))

Parece normal que la mayoría de pisos no tuvieran garaje incluido.

Barrio

Esta variable no será muy útil para el análisis, ya que en muchos casos son descripciones de la vivienda hechas por el usuario, por lo que habría que hacer un tratamiento previo para poder agruparlas. Algo innecesario ya que ese mismo efecto se consigue con el siguiente campo: distritos.

merge(setNames(as.data.frame(table(mhp_train$neighborhood)), c("neighborhood", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$neighborhood))*100, 2)), c("neighborhood", "prop (%)"))
) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
neighborhood count prop (%)
12 de Octubre-Orcasur 19 0.25
Abrantes 54 0.70
Acacias 50 0.65
Adelfas 30 0.39
Aeropuerto 2 0.03
Águilas 33 0.43
Alameda de Osuna 32 0.41
Almagro 140 1.81
Almendrales 42 0.54
Aluche 56 0.72
Ambroz 30 0.39
Amposta 18 0.23
Apóstol Santiago 19 0.25
Arapiles 67 0.87
Aravaca 47 0.61
Arcos 26 0.34
Argüelles 142 1.84
Arroyo del Fresno 11 0.14
Atalaya 6 0.08
Ático en Almagro 8 0.10
Ático en Arapiles 1 0.01
Ático en Aravaca 0 0.00
Ático en Arcos 1 0.01
Ático en Argüelles 0 0.00
Ático en Bernabéu-Hispanoamérica 1 0.01
Ático en Campo de las Naciones-Corralejos 1 0.01
Ático en Castellana 10 0.13
Ático en Castilla 4 0.05
Ático en Chueca-Justicia 1 0.01
Ático en Ciudad Universitaria 1 0.01
Ático en Colina 1 0.01
Ático en Concepción 1 0.01
Ático en Conde Orgaz-Piovera 1 0.01
Ático en Costillares 1 0.01
Ático en Cuatro Caminos 0 0.00
Ático en Delicias 1 0.01
Ático en El Cañaveral 1 0.01
Ático en El Viso 8 0.10
Ático en Fuente del Berro 0 0.00
Ático en Gaztambide 2 0.03
Ático en Goya 3 0.04
Ático en Guindalera 2 0.03
Ático en Horcajo 1 0.01
Ático en Ibiza 2 0.03
Ático en Jerónimos 1 0.01
Ático en Las Tablas 1 0.01
Ático en Lista 1 0.01
Ático en Malasaña-Universidad 1 0.01
Ático en Niño Jesús 1 0.01
Ático en Nueva España 4 0.05
Ático en Nuevos Ministerios-Ríos Rosas 3 0.04
Ático en Pacífico 1 0.01
Ático en Palacio 0 0.00
Ático en Palomas 1 0.01
Ático en Peñagrande 1 0.01
Ático en Pinar del Rey 0 0.00
Ático en Pueblo Nuevo 1 0.01
Ático en Puerta Bonita 0 0.00
Ático en Puerta del Ángel 1 0.01
Ático en Recoletos 3 0.04
Ático en Rejas 0 0.00
Ático en Rosas 1 0.01
Ático en Salvador 1 0.01
Ático en San Juan Bautista 1 0.01
Ático en San Pascual 1 0.01
Ático en Simancas 0 0.00
Ático en Sol 3 0.04
Ático en Trafalgar 1 0.01
Ático en Valdeacederas 1 0.01
Ático en Valdebebas - Valdefuentes 1 0.01
Ático en Valdebernardo - Valderrivas 0 0.00
Ático en Valdemarín 2 0.03
Ático en Vallehermoso 2 0.03
Bellas Vistas 89 1.15
Bernabéu-Hispanoamérica 91 1.18
Berruguete 82 1.06
Buena Vista 40 0.52
Butarque 14 0.18
Campamento 16 0.21
Campo de las Naciones-Corralejos 16 0.21
Canillas 43 0.56
Canillejas 78 1.01
Casa de Campo 18 0.23
Casa o chalet independiente en Aravaca 7 0.09
Casa o chalet independiente en Bernabéu-Hispanoamérica 2 0.03
Casa o chalet independiente en Canillas 2 0.03
Casa o chalet independiente en Canillejas 1 0.01
Casa o chalet independiente en Casa de Campo 0 0.00
Casa o chalet independiente en Casco Histórico de Barajas 0 0.00
Casa o chalet independiente en Castilla 0 0.00
Casa o chalet independiente en Ciudad Universitaria 5 0.06
Casa o chalet independiente en Conde Orgaz-Piovera 7 0.09
Casa o chalet independiente en El Plantío 1 0.01
Casa o chalet independiente en El Viso 4 0.05
Casa o chalet independiente en Fuentelarreina 0 0.00
Casa o chalet independiente en Guindalera 0 0.00
Casa o chalet independiente en Mirasierra 4 0.05
Casa o chalet independiente en Niño Jesús 1 0.01
Casa o chalet independiente en Nueva España 6 0.08
Casa o chalet independiente en Opañel 0 0.00
Casa o chalet independiente en Orcasitas 0 0.00
Casa o chalet independiente en Peñagrande 4 0.05
Casa o chalet independiente en Rejas 1 0.01
Casa o chalet independiente en Tres Olivos - Valverde 2 0.03
Casa o chalet independiente en Valdemarín 7 0.09
Casa o chalet independiente en Valdezarza 0 0.00
Casa o chalet independiente en Villaverde Alto 0 0.00
Casa o chalet independiente en Vista Alegre 0 0.00
Casco Histórico de Barajas 22 0.28
Casco Histórico de Vallecas 87 1.13
Casco Histórico de Vicálvaro 35 0.45
Castellana 172 2.22
Castilla 62 0.80
Chalet adosado en Alameda de Osuna 1 0.01
Chalet adosado en Aravaca 0 0.00
Chalet adosado en Arroyo del Fresno 1 0.01
Chalet adosado en Bernabéu-Hispanoamérica 1 0.01
Chalet adosado en Buena Vista 1 0.01
Chalet adosado en Campo de las Naciones-Corralejos 1 0.01
Chalet adosado en Canillas 1 0.01
Chalet adosado en Canillejas 0 0.00
Chalet adosado en Ciudad Jardín 1 0.01
Chalet adosado en Ciudad Universitaria 0 0.00
Chalet adosado en Concepción 1 0.01
Chalet adosado en Conde Orgaz-Piovera 5 0.06
Chalet adosado en El Plantío 2 0.03
Chalet adosado en El Viso 4 0.05
Chalet adosado en Fuente del Berro 1 0.01
Chalet adosado en Guindalera 0 0.00
Chalet adosado en Los Cármenes 0 0.00
Chalet adosado en Mirasierra 1 0.01
Chalet adosado en Moscardó 1 0.01
Chalet adosado en Nueva España 3 0.04
Chalet adosado en Numancia 0 0.00
Chalet adosado en Orcasitas 2 0.03
Chalet adosado en Palomas 1 0.01
Chalet adosado en Peñagrande 1 0.01
Chalet adosado en Pinar del Rey 1 0.01
Chalet adosado en Puerta Bonita 0 0.00
Chalet adosado en Rejas 1 0.01
Chalet adosado en Salvador 0 0.00
Chalet adosado en Tres Olivos - Valverde 0 0.00
Chalet adosado en Valdeacederas 0 0.00
Chalet adosado en Valdebebas - Valdefuentes 1 0.01
Chalet adosado en Valdezarza 1 0.01
Chalet en Águilas 0 0.00
Chalet en Aravaca 1 0.01
Chalet en Ciudad Universitaria 1 0.01
Chalet en El Plantío 1 0.01
Chalet en El Viso 0 0.00
Chalet en Entrevías 1 0.01
Chalet en Fuentelarreina 0 0.00
Chalet en Mirasierra 1 0.01
Chalet en Rejas 0 0.00
Chalet en Valdemarín 0 0.00
Chalet pareado en Abrantes 1 0.01
Chalet pareado en Alameda de Osuna 0 0.00
Chalet pareado en Apóstol Santiago 1 0.01
Chalet pareado en Aravaca 3 0.04
Chalet pareado en Arroyo del Fresno 1 0.01
Chalet pareado en Bernabéu-Hispanoamérica 0 0.00
Chalet pareado en Campo de las Naciones-Corralejos 0 0.00
Chalet pareado en Canillas 0 0.00
Chalet pareado en Canillejas 0 0.00
Chalet pareado en Ciudad Jardín 0 0.00
Chalet pareado en Ciudad Universitaria 1 0.01
Chalet pareado en Colina 1 0.01
Chalet pareado en Conde Orgaz-Piovera 4 0.05
Chalet pareado en El Plantío 2 0.03
Chalet pareado en El Viso 3 0.04
Chalet pareado en Fuentelarreina 0 0.00
Chalet pareado en Guindalera 0 0.00
Chalet pareado en Mirasierra 3 0.04
Chalet pareado en Nueva España 0 0.00
Chalet pareado en Palomas 4 0.05
Chalet pareado en Peñagrande 1 0.01
Chalet pareado en Rejas 1 0.01
Chalet pareado en Tres Olivos - Valverde 0 0.00
Chalet pareado en Valdebebas - Valdefuentes 1 0.01
Chopera 34 0.44
Chueca-Justicia 67 0.87
Ciudad Jardín 47 0.61
Ciudad Universitaria 40 0.52
Colina 25 0.32
Comillas 52 0.67
Concepción 49 0.63
Conde Orgaz-Piovera 35 0.45
Costillares 53 0.69
Cuatro Caminos 98 1.27
Cuatro Vientos 1 0.01
Cuzco-Castillejos 98 1.27
Delicias 58 0.75
Dúplex en Almagro 6 0.08
Dúplex en Apóstol Santiago 0 0.00
Dúplex en Arapiles 0 0.00
Dúplex en Aravaca 3 0.04
Dúplex en Arcos 0 0.00
Dúplex en Argüelles 0 0.00
Dúplex en Arroyo del Fresno 1 0.01
Dúplex en Bernabéu-Hispanoamérica 0 0.00
Dúplex en Berruguete 1 0.01
Dúplex en Butarque 0 0.00
Dúplex en Canillas 1 0.01
Dúplex en Canillejas 1 0.01
Dúplex en Casco Histórico de Barajas 0 0.00
Dúplex en Casco Histórico de Vallecas 0 0.00
Dúplex en Castellana 1 0.01
Dúplex en Castilla 1 0.01
Dúplex en Chueca-Justicia 1 0.01
Dúplex en Colina 1 0.01
Dúplex en Concepción 2 0.03
Dúplex en Conde Orgaz-Piovera 1 0.01
Dúplex en Costillares 0 0.00
Dúplex en Cuzco-Castillejos 1 0.01
Dúplex en El Viso 4 0.05
Dúplex en Entrevías 1 0.01
Dúplex en Fontarrón 1 0.01
Dúplex en Fuente del Berro 0 0.00
Dúplex en Goya 1 0.01
Dúplex en Imperial 1 0.01
Dúplex en Jerónimos 4 0.05
Dúplex en Las Tablas 2 0.03
Dúplex en Legazpi 2 0.03
Dúplex en Lista 0 0.00
Dúplex en Malasaña-Universidad 0 0.00
Dúplex en Mirasierra 1 0.01
Dúplex en Moscardó 1 0.01
Dúplex en Nueva España 8 0.10
Dúplex en Nuevos Ministerios-Ríos Rosas 1 0.01
Dúplex en Numancia 2 0.03
Dúplex en Opañel 1 0.01
Dúplex en Pacífico 3 0.04
Dúplex en Palacio 1 0.01
Dúplex en Pinar del Rey 0 0.00
Dúplex en Prosperidad 1 0.01
Dúplex en Puerta Bonita 0 0.00
Dúplex en Puerta del Ángel 0 0.00
Dúplex en Quintana 0 0.00
Dúplex en Recoletos 3 0.04
Dúplex en Rejas 0 0.00
Dúplex en San Isidro 1 0.01
Dúplex en San Juan Bautista 6 0.08
Dúplex en Sanchinarro 1 0.01
Dúplex en Simancas 2 0.03
Dúplex en Trafalgar 2 0.03
Dúplex en Tres Olivos - Valverde 1 0.01
Dúplex en Valdeacederas 3 0.04
Dúplex en Valdebernardo - Valderrivas 0 0.00
Dúplex en Valdemarín 4 0.05
Dúplex en Vallehermoso 1 0.01
Dúplex en Ventilla-Almenara 0 0.00
Dúplex en Villaverde Alto 1 0.01
El Cañaveral 50 0.65
El Pardo 3 0.04
El Plantío 3 0.04
El Viso 109 1.41
Ensanche de Vallecas - La Gavia 60 0.78
Entrevías 38 0.49
Estrella 27 0.35
Fontarrón 31 0.40
Fuente del Berro 45 0.58
Fuentelarreina 9 0.12
Gaztambide 77 1.00
Goya 257 3.32
Guindalera 97 1.25
Hellín 17 0.22
Horcajo 3 0.04
Huertas-Cortes 48 0.62
Ibiza 56 0.72
Imperial 48 0.62
Jerónimos 62 0.80
La Paz 48 0.62
Las Tablas 39 0.50
Lavapiés-Embajadores 115 1.49
Legazpi 38 0.49
Lista 134 1.73
Los Ángeles 43 0.56
Los Berrocales 3 0.04
Los Cármenes 50 0.65
Los Rosales 34 0.44
Lucero 43 0.56
Malasaña-Universidad 126 1.63
Marroquina 17 0.22
Media Legua 20 0.26
Mirasierra 15 0.19
Montecarmelo 14 0.18
Moscardó 47 0.61
Niño Jesús 24 0.31
Nueva España 76 0.98
Nuevos Ministerios-Ríos Rosas 91 1.18
Numancia 90 1.16
Opañel 54 0.70
Orcasitas 38 0.49
Pacífico 83 1.07
Palacio 88 1.14
Palomas 20 0.26
Palomeras Bajas 42 0.54
Palomeras sureste 37 0.48
Palos de Moguer 73 0.94
Pau de Carabanchel 26 0.34
Pavones 6 0.08
Peñagrande 77 1.00
Pilar 63 0.81
Pinar del Rey 67 0.87
Portazgo 43 0.56
Pradolongo 32 0.41
Prosperidad 86 1.11
Pueblo Nuevo 132 1.71
Puerta Bonita 89 1.15
Puerta del Ángel 137 1.77
Quintana 68 0.88
Recoletos 181 2.34
Rejas 59 0.76
Rosas 30 0.39
Salvador 21 0.27
San Cristóbal 30 0.39
San Diego 120 1.55
San Fermín 33 0.43
San Isidro 91 1.18
San Juan Bautista 39 0.50
San Pascual 48 0.62
Sanchinarro 61 0.79
Santa Eugenia 16 0.21
Simancas 60 0.78
Sol 66 0.85
Timón 14 0.18
Trafalgar 97 1.25
Tres Olivos - Valverde 51 0.66
Valdeacederas 104 1.35
Valdebebas - Valdefuentes 42 0.54
Valdebernardo - Valderrivas 22 0.28
Valdemarín 19 0.25
Valdezarza 42 0.54
Vallehermoso 72 0.93
Ventas 103 1.33
Ventilla-Almenara 46 0.59
Villaverde Alto 87 1.13
Vinateros 20 0.26
Virgen del Cortijo - Manoteras 31 0.40
Vista Alegre 102 1.32
Zofío 25 0.32

Distrito

En este caso los datos si están correctamente recogidos, repartiendo las casas entre los 21 distritos de Madrid.Casi todos con una representación considerable, en lo que destaca el barrio de Salamanca, con casi el doble de viviendas que el segundo que más tiene, Chamberí.

merge(setNames(as.data.frame(table(mhp_train$district)), c("district", "count")),
      setNames(as.data.frame(round(prop.table(table(mhp_train$district))*100, 2)), c("district", "prop (%)"))
) %>%
  arrange(desc(count)) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
district count prop (%)
barrio de salamanca 911 11.78
chamberi 571 7.38
ciudad lineal 540 6.98
chamartin 526 6.80
tetuan 523 6.76
centro 517 6.69
carabanchel 512 6.62
puente-de-vallecas 374 4.84
fuencarral 356 4.60
moncloa 353 4.57
hortaleza 352 4.55
latina 337 4.36
san-blas 319 4.13
Arganzuela 305 3.94
retiro 295 3.82
usera 240 3.10
villaverde 209 2.70
villa-de-vallecas 163 2.11
vicalvaro 141 1.82
moratalaz 99 1.28
barajas 89 1.15
ggplot(mhp_train, aes(district)) +
  geom_bar(fill = "#0BB363") +
  labs(x = "Distrito", y = "Número de viviendas", title = "Viviendas por distrito") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

3.2 Análisis de variables cuantitativas

Precio

data.frame(summarise(mhp_train,
                     min = min(price),
                     max = max(price),
                     median = median(price),
                     mean = mean(price),
                     sd = sd(price))) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
min max median mean sd
47500 8950000 369000 631920.1 753136.6
ggplot(mhp_train, aes(y = price)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Precio", title = "Boxplot de precios") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(price)) +
  geom_histogram(aes(y = ..count..), bins = 50, position = "dodge", fill = "#0BB363") +
  labs(x = "Precio", y = "Número de viviendas", title = "Viviendas por Precio") +
  scale_x_continuous(breaks = pretty(mhp_train$price, n = 10), labels = comma_format()) +
  scale_y_continuous(labels = function(x) sprintf("%.0f", x)) +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30))

Habitaciones

data.frame(summarise(mhp_train,
                     min = min(rooms),
                     max = max(rooms),
                     median = median(rooms),
                     mean = mean(rooms),
                     sd = sd(rooms))) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
min max median mean sd
1 17 3 2.848552 1.302995
ggplot(mhp_train, aes(y = rooms)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Habitaciones", title = "Boxplot de habitaciones") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(rooms)) +
  geom_histogram(aes(), bins = 16, position = "dodge", fill = "#0BB363") +
  geom_density(alpha=.2, fill = "red") +
  labs(x = "Habitaciones", y = "Número de viviendas", title = "Viviendas por número de habitaciones") +
  theme(plot.title = element_text(hjust = 0.5))

data.frame(summarise(mhp_train,
                     min = min(m2),
                     max = max(m2),
                     median = median(m2),
                     mean = mean(m2),
                     sd = sd(m2))) %>%
  kbl() %>%
  kable_material(c("striped", "hover"))
min max median mean sd
20 986 95 127.1208 102.497
ggplot(mhp_train, aes(y = m2)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Metros cuadrados", title = "Boxplot de superficies") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(m2)) +
  geom_histogram(aes(y=..count..), bins = 50, position = "dodge", fill = "#0BB363") +
  geom_density(alpha=.3, fill = "red") +
  labs(x = "Metros cuadrados", y = "Número de viviendas", title = "Viviendas por metros cuadrados") +
  theme(plot.title = element_text(hjust = 0.5))

4 Análisis multivariante

Una vez analizadas todas las variables de forma individual, podemos buscar algunas relaciones entre ellas.

Como se ve en los gráficos inferiores hay una relción fuerte entre el precio, los metros cuadradados y el número de habitaciones.

plot_grid(
  ggcorrplot(cor(mhp_train %>% select_if(is.numeric)),type = "lower", lab=TRUE),
  
  ggplot(mhp_train, aes(x = m2, y = price)) +
  geom_point() +
  geom_smooth() +
  ggtitle('Reación precio y m²') +
  theme(plot.title = element_text(hjust = 0.5)),
  
  ggplot(mhp_train, aes(x = rooms, y = price)) +
  geom_point() +
  geom_smooth() +
  ggtitle('Reación precio y habitaciones') +
  theme(plot.title = element_text(hjust = 0.5)), 
  
  ggplot(mhp_train, aes(x = rooms, y = m2)) +
  geom_point() +
  geom_smooth() +
  ggtitle('Reación m² y habitaciones') +
  theme(plot.title = element_text(hjust = 0.5)), 
  
  nrow = 2
)

Introduciendo algunas de las variables categóricas se puede observar como varía el precio con la superficie o el número de habitaciones dependiendo de si incluye garaje, ascensor u orientación exterior.

plot_grid(
  ggplot(mhp_train, aes(x = m2, y = price, colour = exterior)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-superficie, exterior') ,
  
  ggplot(mhp_train, aes(x = rooms, y = price, colour = exterior)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) + 
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-habitaciones, exterior') ,
  
  nrow =1
)

plot_grid(
  ggplot(mhp_train, aes(x = m2, y = price, colour = elevator)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-superficie, ascensor'),
  
  ggplot(mhp_train, aes(x = rooms, y = price, colour = elevator)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-habitación, ascensor'),
  
  nrow =1
)

plot_grid(
  ggplot(mhp_train, aes(x = m2, y = price, colour = garage)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-superficie, garaje'),
  
  ggplot(mhp_train, aes(x = rooms, y = price, colour = garage)) +
  geom_point() +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  ggtitle('Rel. precio-habitaciones, garaje'),
  
  nrow =1
)

Representando boxplot por distrito se ve claramente que los barrios ricos como Salamanca o Chamberí, a parte de tener precios medios más altos, tienen mayor número de viviendas con precios atípicos (por encima), mientras que en el número de habitaciones no se aprecia una diferencia tan notoria.

ggplot(mhp_train, aes(district, price)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Precio", title = "Boxplot de precios por distrito") +
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(district, rooms)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Habitaciones", title = "Boxplot de habitaciones por distrito") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

ggplot(mhp_train, aes(district, m2)) +
  geom_boxplot(fill = "#0BB363") +
  labs(y = "Metros cuadrados", title = "Boxplot de superficie por distrito") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1))

En las siguientes figuras se representa el precio medio por distrito y como varía con respecto a otras variables de interés.

mhp_train %>%
  group_by(district, m2) %>%
  summarize(avg_price = mean(price)) %>%
  ggplot(aes(x = m2, y = avg_price)) + 
  geom_point(size = 0.5) +
  facet_wrap(~ district) + 
  scale_y_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6))

mhp_train %>%
  group_by(district, exterior) %>% 
  summarise(avg_price = mean(price)) %>%
  ggplot(aes(x=district, y=avg_price, fill=exterior)) +
  geom_bar(stat = "identity", position = "dodge") + 
  scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  ggtitle("Precio medio por distrito y exterior") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) + 
  labs(x = "Distrito", y = "Precio medio") + 
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))

mhp_train %>%
  group_by(district, elevator) %>% 
  summarise(avg_price = mean(price)) %>%
  ggplot(aes(x=district, y=avg_price, fill=elevator)) +
  geom_bar(stat = "identity", position = "dodge") + 
  scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  ggtitle("Precio medio por distrito y ascensor") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) + 
  labs(x = "Distrito", y = "Precio medio") + 
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))

mhp_train %>%
  group_by(district, garage) %>% 
  summarise(avg_price = mean(price)) %>%
  ggplot(aes(x=district, y=avg_price, fill=garage)) +
  geom_bar(stat = "identity", position = "dodge") + 
  scale_fill_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  ggtitle("Precio medio por distrito y garaje") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 30, hjust = 1)) + 
  labs(x = "Distrito", y = "Precio medio") + 
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))

5 Procesado de variables cualitativas

mhp_train %>%
  group_by(house_type) %>%
  summarise(precio_medio = mean(price)) %>%
  arrange(precio_medio)
## # A tibble: 26 × 2
##    house_type  precio_medio
##    <fct>              <dbl>
##  1 sotano           279031.
##  2 planta -1        298188.
##  3 semi-sotano      302916.
##  4 bajo             316306.
##  5 planta 11        458019.
##  6 planta 1         540375.
##  7 planta 2         595355.
##  8 entreplanta      598834.
##  9 planta 10        605368.
## 10 planta 4         605711.
## # … with 16 more rows

Para el posterior análisis transformamos la columna house_type, agrupando en “piso” (independientemente de la planta en la que esté), “casa” (casa o chalet) y “sotano”.

mhp_train <- mhp_train %>%
  mutate(tipo_casa = case_when(
    grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
    grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
    TRUE ~ "piso"
  ))
mhp_train$tipo_casa <- as.factor(mhp_train$tipo_casa)

También se agruparán los distritos en función del precio medio del mismo: caro, medio o barato.

mhp_train <- mhp_train %>%
  group_by(district) %>%
  mutate(precio_distrito = mean(price)) %>%
  ungroup()

# Calcular los percentiles
percentiles <- quantile(mhp_train$precio_distrito, probs = c(0, 1/3, 2/3, 1))

# Asignar etiquetas
mhp_train$precio_distrito <- cut(mhp_train$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)

Además, de cara a desarrollar algunos modelos, se creará una variable target binaria que dividirá las casas entre “cara” y “barata” (1/0).

mhp_train$precio_bin <- ifelse(mhp_train$price > median(mhp_train$price), "cara", "barata")
mhp_train$precio_bin = ifelse(mhp_train$precio_bin == "cara", 1, 0)
mhp_train$precio_bin = factor(mhp_train$precio_bin, levels = c(1,0))
# Aplicamos todos los cambios también a test.

# Categorizamos house_type
mhp_test <- mhp_test %>%
  mutate(tipo_casa = case_when(
    grepl("casa|chalet", house_type, ignore.case = TRUE) ~ "casa",
    grepl("\\bplanta\\s*-?1\\b|semi-?sotano|bajo|sotano", house_type, ignore.case = TRUE) ~ "sotano",
    TRUE ~ "piso"
  ))
mhp_test$tipo_casa <- as.factor(mhp_test$tipo_casa)

# Categorizamos los distritos
mhp_test <- mhp_test %>%
  group_by(district) %>%
  mutate(precio_distrito = mean(price)) %>%
  ungroup()

percentiles <- quantile(mhp_test$precio_distrito, probs = c(0, 1/3, 2/3, 1))

mhp_test$precio_distrito <- cut(mhp_test$precio_distrito, breaks = percentiles, labels = c("barato", "medio", "caro"), include.lowest = TRUE)

# Creamos variable target
mhp_test$precio_bin <- ifelse(mhp_test$price > median(mhp_test$price), "cara", "barata")
mhp_test$precio_bin = ifelse(mhp_test$precio_bin == "cara", 1, 0)
mhp_test$precio_bin = factor(mhp_test$precio_bin, levels = c(1,0))

6 Modelo de regresión lineal

Con el dataset ya preparado, se crea un modelo de regresión lineal.

En este caso, tras hacer diversas pruebas se ha optado por introducir sólo 3 variables (tipo_casa, m2 y precio_distrito), alcanzando un R2 de 0.7553.

Esta cifra podía incrementarse algo, sin embargo es a costa de introducir más variables y por tanto más complejidad al modelo, haciéndolo menos explicable.

lm_fit <- lm(price ~ m2+tipo_casa+precio_distrito, data=mhp_train)
summary(lm_fit)
## 
## Call:
## lm(formula = price ~ m2 + tipo_casa + precio_distrito, data = mhp_train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3277290  -137875   -17495    95705  5483233 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -915949.86   30574.68 -29.958   <2e-16 ***
## m2                      6052.79      51.28 118.027   <2e-16 ***
## tipo_casapiso         673879.37   27865.36  24.183   <2e-16 ***
## tipo_casasotano       606283.79   28773.95  21.071   <2e-16 ***
## precio_distritomedio  101183.27   10337.15   9.788   <2e-16 ***
## precio_distritocaro   373747.30   11419.31  32.729   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 372600 on 7726 degrees of freedom
## Multiple R-squared:  0.7553, Adjusted R-squared:  0.7552 
## F-statistic:  4771 on 5 and 7726 DF,  p-value: < 2.2e-16
coef(lm_fit)
##          (Intercept)                   m2        tipo_casapiso 
##          -915949.864             6052.795           673879.371 
##      tipo_casasotano precio_distritomedio  precio_distritocaro 
##           606283.787           101183.273           373747.298
residuals=lm_fit$residuals
autoplot(lm_fit)

Observando los residuos se ve como hay cierta heterocedasticidad y en la gráfica QQ no se ajustan a la normal, por lo que el módelo no es bueno.

# Realizar predicciones en la partición test
predictions <- predict(lm_fit, newdata = mhp_test)

# Comparar las predicciones con los valores reales
rmse <- rmse(predictions, mhp_test$price)
mae <- mae(predictions, mhp_test$price)
r2 <- R2(predictions, mhp_test$price)

# Imprimir las métricas de evaluación
cat(paste0("RMSE: ", round(rmse, 2), "\n"))
## RMSE: 398263.79
cat(paste0("MAE: ", round(mae, 2), "\n"))
## MAE: 199633.64
cat(paste0("R-squared: ", round(r2, 2)))
## R-squared: 0.73

Un R2 de 0.73 (bastante similar al calculado en train) indica que se explica con esa probabilidad la variabilidad de la variable dependiente, por lo tanto podríamos considerarlo razonablemente útil, sin embargo sería preferible mejorar los valores del error de predicción promedio o el cuadrático.

7 Conclusiones preliminares

El análisis realizado hasta el momento ha servido para familiarizarse con el mercado inmobiliario de Madrid (o al menos parte de él) y conocer como se comportan algunas de las variables más relevantes en el valor de las viviendas.

Con estos conocimientos se abordarán nuevos modelos para intentar conseguir mejores predicciones y simultaneamente comprender y desarrollar estos mismos modelos.

str(mhp_train)
## tibble [7,732 × 12] (S3: tbl_df/tbl/data.frame)
##  $ price          : num [1:7732] 499000 385000 950000 460000 335000 580000 940000 85300 399000 125000 ...
##  $ house_type     : Factor w/ 26 levels "bajo","casa",..: 21 16 20 10 18 19 16 16 6 1 ...
##  $ exterior       : Factor w/ 2 levels "1","0": 1 1 1 1 1 1 1 1 1 1 ...
##  $ rooms          : num [1:7732] 4 1 3 3 4 3 2 3 3 2 ...
##  $ m2             : num [1:7732] 131 60 160 114 122 120 124 75 86 62 ...
##  $ elevator       : Factor w/ 2 levels "1","0": 1 2 1 1 1 1 1 1 2 2 ...
##  $ garage         : Factor w/ 2 levels "1","0": 2 2 2 1 2 2 2 2 2 2 ...
##  $ neighborhood   : Factor w/ 341 levels "12 de Octubre-Orcasur",..: 76 17 8 284 318 265 264 316 192 11 ...
##  $ district       : Factor w/ 21 levels "Arganzuela","barajas",..: 17 12 7 13 18 3 3 21 1 19 ...
##  $ tipo_casa      : Factor w/ 3 levels "casa","piso",..: 2 2 2 2 2 2 2 2 3 3 ...
##  $ precio_distrito: Factor w/ 3 levels "barato","medio",..: 2 3 3 1 1 3 3 1 1 1 ...
##  $ precio_bin     : Factor w/ 2 levels "1","0": 1 1 1 1 2 1 1 2 1 2 ...

8 Aprendizaje no supervisado

8.1 K-Means

El algoritmo K-Means es un método de aprendizaje automático no supervisado utilizado para el análisis de agrupamiento o “clustering”. El objetivo principal de este algoritmo es dividir los datos en K grupos (clusters), de modo que los puntos dentro de un mismo grupo sean lo más similares posible.

Por como funciona este método, es habitual utilizarlo en la etapa de EDA para agrupar y etiquetar observaciones con un criterior más científico que el basado en el conocimiento personal del dominio, más que para hacer predicciones (aunque también puede hacerse).

En nuestro lo utilizaremos para comprobar si hay alguna forma mejor de agrupar las casas de como lo hemos hecho anteriormente.

mhp_train_kmeans <- subset(mhp_train, select = c("m2", "rooms", "elevator", "garage", "exterior"))

mhp_train_kmeans$exterior <- as.numeric(mhp_train_kmeans$exterior)
mhp_train_kmeans$elevator <- as.numeric(mhp_train_kmeans$elevator)
mhp_train_kmeans$garage <- as.numeric(mhp_train_kmeans$garage)

# Normalizamos los datos
mhp_train_kmeans <- scale(mhp_train_kmeans)

Tras normalizar los datos, aplicamos k-means sólo para las variables numéricas y mediante los siguientes test decidiremos cual es el número óptimo de clusters.

fviz_nbclust(mhp_train_kmeans, kmeans, method = "wss")

fviz_nbclust(mhp_train_kmeans, kmeans, method = "silhouette")

Utilizando la regla del codo en el gráfico wss podríamos decir que está en el 6, algo que confirma el segundo gráfico. Por lo tanto, aplicamos el modelo con K = 6

# Definir el número de grupos
k <- 6

# Ejecutar k-means
grupos <- kmeans(mhp_train_kmeans, k)

mhp_train$cluster <- as.factor(grupos$cluster)

# Graficar la relación entre las características y los grupos
ggplot(mhp_train, aes(x=m2, y=price, color=cluster)) + 
  geom_point() +
  ggtitle("Grupos de casas en venta") +
  theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Tamaño en m2") +
  ylab("Precio")

fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "euclid",repel = TRUE,star.plot = TRUE)

fviz_cluster(grupos, data = mhp_train_kmeans, ellipse.type = "norm")

Graficamente es difícil diferenciar las agrupaciones, por lo que en la tabla inferior obtenemos los precios medios de cada cluster, donde se puede apreciar que la gran mayoría de agrupaciones son coherentes (salvo tal vez el 2 y el 3 que son muy similares).

También se muestra el tipo de casa y el distrito predominante en cada cluster.

# Precio medio por cluster
mhp_train %>%
  group_by(cluster) %>%
  summarise(precio_medio = mean(price))
## # A tibble: 6 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            184203.
## 2 2            640443.
## 3 3            617779.
## 4 4            223476.
## 5 5           2228444.
## 6 6            355840.
# Tipo de casa predominante por cluster
mhp_train %>%
  group_by(cluster, tipo_casa) %>%
  summarise(frecuencia = n()) %>%
  arrange(cluster, desc(frecuencia)) %>%
  slice(1) %>%
  select(cluster, tipo_casa_predominante = tipo_casa)
## # A tibble: 6 × 2
## # Groups:   cluster [6]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       sotano                
## 2 2       piso                  
## 3 3       piso                  
## 4 4       piso                  
## 5 5       piso                  
## 6 6       piso
# Distrito predominante por cluster
mhp_train %>%
  group_by(cluster, district) %>%
  summarise(frecuencia = n()) %>%
  arrange(cluster, desc(frecuencia)) %>%
  slice(1) %>%
  select(cluster, distrito_predominante = district)
## # A tibble: 6 × 2
## # Groups:   cluster [6]
##   cluster distrito_predominante
##   <fct>   <fct>                
## 1 1       tetuan               
## 2 2       hortaleza            
## 3 3       barrio de salamanca  
## 4 4       carabanchel          
## 5 5       barrio de salamanca  
## 6 6       barrio de salamanca

Sabemos que tanto la regla del codo como cualquier otro método para definir el número de clusters no siempre es adecuada, por lo que a continuación se prueba con un bucle que prueba con número de clusters menores que 6 y muestra las agrupaciones por precio y tipo de piso anteriores.

# Ejecuta k-means para diferentes valores de K
for (k in 2:6) {
  # Ejecutar k-means
  grupos <- kmeans(mhp_train_kmeans, centers = k)
  mhp_train$cluster <- as.factor(grupos$cluster)
  
  print(k)
  
  print(mhp_train %>%
  group_by(cluster) %>%
  summarise(precio_medio = mean(price)))
  
  print(mhp_train %>%
  group_by(cluster, tipo_casa) %>%
  summarise(frecuencia = n()) %>%
  arrange(cluster, desc(frecuencia)) %>%
  slice(1) %>%
  select(cluster, tipo_casa_predominante = tipo_casa))
}
## [1] 2
## # A tibble: 2 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            394442.
## 2 2           1094473.
## # A tibble: 2 × 2
## # Groups:   cluster [2]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## [1] 3
## # A tibble: 3 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            402986.
## 2 2            605000.
## 3 3           1941976.
## # A tibble: 3 × 2
## # Groups:   cluster [3]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       piso                  
## [1] 4
## # A tibble: 4 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            579636.
## 2 2            352563.
## 3 3            665255.
## 4 4            189057.
## # A tibble: 4 × 2
## # Groups:   cluster [4]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       piso                  
## 4 4       sotano                
## [1] 5
## # A tibble: 5 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1            435798.
## 2 2            241929.
## 3 3            189057.
## 4 4            665853.
## 5 5           1005091.
## # A tibble: 5 × 2
## # Groups:   cluster [5]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       sotano                
## 4 4       piso                  
## 5 5       piso                  
## [1] 6
## # A tibble: 6 × 2
##   cluster precio_medio
##   <fct>          <dbl>
## 1 1           1179261.
## 2 2            652498.
## 3 3            355840.
## 4 4            217255.
## 5 5           2356279.
## 6 6            451222.
## # A tibble: 6 × 2
## # Groups:   cluster [6]
##   cluster tipo_casa_predominante
##   <fct>   <fct>                 
## 1 1       piso                  
## 2 2       piso                  
## 3 3       piso                  
## 4 4       piso                  
## 5 5       piso                  
## 6 6       piso

Observando las diferencias entre los diferentes modelos, el número de clusters que mejor diferencia las casas por su precio es 3, algo que encaja con la presunción que inicial que hicimos para dividir los tipos de casa en 3 tipos.

Si observamos además el tipo de casas más común por cluster, casi siempre es piso (salvo en los casos de clusters con muy poco precio medio, que es sótano), lo cual aunque tiene sentido porque este tipo es claramente mayoritario en todo el dataset, no nos permite apreciar a simple vista donde se acumulan el resto de viviendas.

Por lo tanto, entre los datos extraídos del análisis de clusters y el conocimiento previo que teníamos del dominio, tiene sentido mantener los 3 tipos de viviendas principales como variable para los siguientes modelos (sótano, piso y casa).

9 Técnicas de reducción de la dimensionalidad

9.1 PCA

El Análisis de Componentes Principales es una técnica de reducción de la dimensionalidad que se utiliza para transformar un conjunto de datos de múltiples variables en un nuevo conjunto de datos con un número menor de variables cuyo principal objetivo es simplificar la representación de datos mientras se conserva la mayor cantidad de información y variabilidad posible identificando patrones y estructuras subyacentes, reduciendo la multicolinealidad y el sobreajuste.

PAra nuestro dataset, seleccionaremos las columnas numéricas, descartando la columna precio puesto que será de alguna manera nuestra variable dependiente. Además creamos 3 variables dummy usando variables categóricas de tipo 0|1.

train_pca <- mhp_train[,c(3:7,10,11,12)]
dummy_garaje <- model.matrix(~0 + garage,data = train_pca)
dummy_exterior <- model.matrix(~0 + exterior, data = train_pca)
dummy_ascensor <- model.matrix(~0 + elevator, data = train_pca)
dummy_tipo_casa <- model.matrix(~0 + tipo_casa, data = train_pca)
dummy_distrito <- model.matrix(~0 + precio_distrito, data = train_pca)
dummy_data <- cbind(train_pca[, c("rooms", "m2")],dummy_exterior,dummy_garaje,dummy_ascensor,dummy_distrito,dummy_tipo_casa)
prcomp(dummy_data)
## Standard deviations (1, .., p=14):
##  [1] 1.025015e+02 9.457105e-01 7.267966e-01 6.257613e-01 6.127629e-01
##  [6] 5.782676e-01 4.571490e-01 4.037562e-01 1.716715e-01 2.268386e-15
## [11] 1.853843e-15 1.003401e-15 8.037840e-16 4.438592e-16
## 
## Rotation (n x k) = (14 x 14):
##                                 PC1           PC2          PC3          PC4
## rooms                 -8.797027e-03 -0.9916184565 -0.014616928  0.094846395
## m2                    -9.999557e-01  0.0085416356  0.001971791 -0.001450128
## exterior1             -5.771447e-04 -0.0453751425 -0.009449951 -0.063478191
## exterior0              5.771447e-04  0.0453751425  0.009449951  0.063478191
## garage1               -1.538147e-03  0.0421269404 -0.223561300  0.037355645
## garage0                1.538147e-03 -0.0421269404  0.223561300 -0.037355645
## elevator1             -4.274374e-04  0.0241411854 -0.476833313  0.342257253
## elevator0              4.274374e-04 -0.0241411854  0.476833313 -0.342257253
## precio_distritobarato  1.400700e-03 -0.0385563529  0.200000119 -0.359200946
## precio_distritomedio   6.834450e-05 -0.0051663464 -0.051807433  0.201186914
## precio_distritocaro   -1.469044e-03  0.0437226993 -0.148192687  0.158014031
## tipo_casacasa         -8.085078e-04  0.0005892036  0.067961955 -0.037840260
## tipo_casapiso          2.818685e-05 -0.0469938502 -0.466029098 -0.505128722
## tipo_casasotano        7.803209e-04  0.0464046466  0.398067143  0.542968982
##                                 PC5           PC6          PC7          PC8
## rooms                 -0.0132115607 -0.0291253564  0.074126853 -0.029792435
## m2                     0.0009596296  0.0005665667 -0.002313108 -0.001386419
## exterior1             -0.1552630611 -0.0333936354 -0.327138796  0.601052285
## exterior0              0.1552630611  0.0333936354  0.327138796 -0.601052285
## garage1               -0.5627134868 -0.1309300044  0.333444368  0.033280690
## garage0                0.5627134868  0.1309300044 -0.333444368 -0.033280690
## elevator1              0.0482119831 -0.1533456743 -0.318624110 -0.141577832
## elevator0             -0.0482119831  0.1533456743  0.318624110  0.141577832
## precio_distritobarato -0.2366848900 -0.4502445095 -0.365221920 -0.322493702
## precio_distritomedio  -0.1969164022  0.7585991702 -0.088292084 -0.036423761
## precio_distritocaro    0.4336012922 -0.3083546607  0.453514004  0.358917463
## tipo_casacasa         -0.0381883034  0.0229829599  0.061657352 -0.001868362
## tipo_casapiso          0.1388358808  0.1326235059  0.004437444 -0.006913037
## tipo_casasotano       -0.1006475774 -0.1556064657 -0.066094796  0.008781399
##                                PC9          PC10          PC11          PC12
## rooms                  0.000862928  5.877890e-16  1.174221e-15  7.345674e-17
## m2                    -0.001036478 -4.349751e-20  3.574359e-18 -4.762673e-18
## exterior1              0.017792902  8.077407e-03  2.416838e-02 -9.348594e-03
## exterior0             -0.017792902  8.077407e-03  2.416838e-02 -9.348594e-03
## garage1               -0.027676182  6.464949e-03  5.102210e-03 -8.893342e-02
## garage0                0.027676182  6.464949e-03  5.102210e-03 -8.893342e-02
## elevator1              0.086614598  6.353245e-01 -3.028358e-01 -6.708710e-02
## elevator0             -0.086614598  6.353245e-01 -3.028358e-01 -6.708710e-02
## precio_distritobarato -0.004863213 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritomedio  -0.010432322 -9.943406e-02 -8.213108e-02 -5.584679e-01
## precio_distritocaro    0.015295535 -9.943406e-02 -8.213108e-02 -5.584679e-01
## tipo_casacasa          0.809211511 -2.329879e-01 -5.148217e-01  1.145226e-01
## tipo_casapiso         -0.393939361 -2.329879e-01 -5.148217e-01  1.145226e-01
## tipo_casasotano       -0.415272150 -2.329879e-01 -5.148217e-01  1.145226e-01
##                                PC13          PC14
## rooms                 -9.651041e-17 -5.224615e-18
## m2                     4.139246e-18 -4.009145e-19
## exterior1              5.278108e-02 -7.046115e-01
## exterior0              5.278108e-02 -7.046115e-01
## garage1               -6.995888e-01 -5.097578e-02
## garage0               -6.995888e-01 -5.097578e-02
## elevator1              1.228500e-02 -1.293890e-03
## elevator0              1.228500e-02 -1.293890e-03
## precio_distritobarato  6.884856e-02  8.609917e-03
## precio_distritomedio   6.884856e-02  8.609917e-03
## precio_distritocaro    6.884856e-02  8.609917e-03
## tipo_casacasa         -1.877163e-02 -2.325502e-02
## tipo_casapiso         -1.877163e-02 -2.325502e-02
## tipo_casasotano       -1.877163e-02 -2.325502e-02

Podemos observar que en la PC1 lleva un peso casi de 100 la variable m2, mientras que en PC2 pasa lo mismo con la variable rooms. En el resto de PCs el análisis de componentes principales está algo más repartido.

summary(prcomp(dummy_data))
## Importance of components:
##                             PC1     PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     102.5015 0.94571 0.72680 0.62576 0.61276 0.57827 0.45715
## Proportion of Variance   0.9997 0.00009 0.00005 0.00004 0.00004 0.00003 0.00002
## Cumulative Proportion    0.9997 0.99981 0.99986 0.99989 0.99993 0.99996 0.99998
##                            PC8    PC9      PC10      PC11      PC12      PC13
## Standard deviation     0.40376 0.1717 2.268e-15 1.854e-15 1.003e-15 8.038e-16
## Proportion of Variance 0.00002 0.0000 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  1.00000 1.0000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
##                             PC14
## Standard deviation     4.439e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

Al hacer un summary de las variables vemos que la proporción de varianza explicada por la primera componente es de un 99%, es decir, que es la única relevante. Si nos fijamos vemos que la PC1 está explicada casi al 100% por los m2. En nuestra tabla vemos que los m2 tienen valores mucho más altos que el resto de variables. Por lo tanto, lo que hacemos a continuación es un análisis pero escalando las variables.

prcomp(dummy_data, scale = T)
## Standard deviations (1, .., p=14):
##  [1] 1.788828e+00 1.554767e+00 1.391248e+00 1.320487e+00 1.255039e+00
##  [6] 1.218340e+00 9.006738e-01 7.679348e-01 4.930742e-01 2.297198e-14
## [11] 4.169026e-15 2.788355e-15 1.325148e-15 8.186691e-16
## 
## Rotation (n x k) = (14 x 14):
##                               PC1          PC2         PC3         PC4
## rooms                 -0.32118124  0.185174155  0.19374853 -0.18231760
## m2                    -0.39188093  0.178866531  0.28992752 -0.15926158
## exterior1             -0.22525556  0.323113180 -0.41045579  0.16000993
## exterior0              0.22525556 -0.323113180  0.41045579 -0.16000993
## garage1               -0.39123327  0.104563202  0.02079773  0.20853492
## garage0                0.39123327 -0.104563202 -0.02079773 -0.20853492
## elevator1             -0.28418360 -0.450986811 -0.07378472  0.27110372
## elevator0              0.28418360  0.450986811  0.07378472 -0.27110372
## precio_distritobarato  0.20341246  0.208545636 -0.33198104  0.01261772
## precio_distritomedio  -0.02774634 -0.002015776  0.03875669  0.15703618
## precio_distritocaro   -0.18188148 -0.213762808  0.30358073 -0.17520956
## tipo_casacasa         -0.14205765  0.358486679  0.29646454 -0.20249803
## tipo_casapiso         -0.16858177 -0.251362769 -0.39199358 -0.48519863
## tipo_casasotano        0.22431380  0.124024528  0.29030539  0.56945631
##                               PC5          PC6          PC7          PC8
## rooms                  0.17534372 -0.091653082  0.472849806  0.496683083
## m2                     0.12920863 -0.038990789  0.209138336  0.106709918
## exterior1              0.31433452 -0.162497819 -0.145528928 -0.064746601
## exterior0             -0.31433452  0.162497819  0.145528928  0.064746601
## garage1               -0.40131726  0.276379019 -0.203710793  0.095975695
## garage0                0.40131726 -0.276379019  0.203710793 -0.095975695
## elevator1              0.11627307  0.012675714  0.278402861 -0.212854201
## elevator0             -0.11627307 -0.012675714 -0.278402861  0.212854201
## precio_distritobarato -0.04064647  0.512859981  0.434958282 -0.014241486
## precio_distritomedio  -0.38476634 -0.695032433  0.035559829  0.020880654
## precio_distritocaro    0.43936573  0.186858173 -0.486898651 -0.006820754
## tipo_casacasa         -0.08089606 -0.001818982  0.145439408 -0.763400011
## tipo_casapiso         -0.15088701 -0.039033227 -0.059791130  0.108495318
## tipo_casasotano        0.18370301  0.040471673  0.007300261  0.171061449
##                               PC9          PC10          PC11          PC12
## rooms                 -0.53136448  1.871326e-14  6.573067e-16 -8.713797e-16
## m2                     0.79477851  5.562477e-15  5.135097e-16  5.962879e-16
## exterior1             -0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## exterior0              0.01407624 -1.579692e-03 -7.045930e-01 -5.852591e-02
## garage1               -0.06271669 -7.002484e-04 -1.864879e-03  2.225429e-02
## garage0                0.06271669 -7.002484e-04 -1.864879e-03  2.225429e-02
## elevator1             -0.02034701 -1.325824e-03 -1.363473e-02  3.145845e-02
## elevator0              0.02034701 -1.325824e-03 -1.363473e-02  3.145845e-02
## precio_distritobarato  0.09621994 -7.055211e-02  4.768696e-02 -5.771886e-01
## precio_distritomedio   0.02039646 -7.038638e-02  4.757494e-02 -5.758327e-01
## precio_distritocaro   -0.12064810 -6.816645e-02  4.607446e-02 -5.576714e-01
## tipo_casacasa         -0.22215925 -2.483203e-01 -1.912383e-03  3.007261e-02
## tipo_casapiso          0.04299775 -6.861954e-01 -5.284580e-03  8.310109e-02
## tipo_casasotano        0.03813232 -6.729672e-01 -5.182706e-03  8.149911e-02
##                                PC13          PC14
## rooms                 -4.757936e-17  3.988302e-17
## m2                    -1.440589e-18 -3.087024e-16
## exterior1              1.099278e-02  3.536496e-04
## exterior0              1.099278e-02  3.536496e-04
## garage1                2.157713e-02 -7.064242e-01
## garage0                2.157713e-02 -7.064242e-01
## elevator1             -7.059752e-01 -2.053508e-02
## elevator0             -7.059752e-01 -2.053508e-02
## precio_distritobarato -2.595459e-02 -1.903173e-02
## precio_distritomedio  -2.589362e-02 -1.898702e-02
## precio_distritocaro   -2.507696e-02 -1.838819e-02
## tipo_casacasa          1.806856e-03  1.253756e-03
## tipo_casapiso          4.992972e-03  3.464565e-03
## tipo_casasotano        4.896719e-03  3.397776e-03

Tras escalar las variables podemos observar que están mucho más repartidos los pesos en las componentes principales debido a que estás asegurando que todas las variables tengan la misma importancia en el análisis, independientemente de sus escalas originales y unidades de medida. Con esto evitas el dominio de variables con mayor magnitud y mejoras la interpretación de resultados.

pca_result <- prcomp(dummy_data, center = TRUE, scale = TRUE)
summary(pca_result)
## Importance of components:
##                           PC1    PC2    PC3    PC4    PC5    PC6     PC7
## Standard deviation     1.7888 1.5548 1.3912 1.3205 1.2550 1.2183 0.90067
## Proportion of Variance 0.2286 0.1727 0.1383 0.1245 0.1125 0.1060 0.05794
## Cumulative Proportion  0.2286 0.4012 0.5395 0.6640 0.7765 0.8826 0.94051
##                            PC8     PC9      PC10      PC11      PC12      PC13
## Standard deviation     0.76793 0.49307 2.297e-14 4.169e-15 2.788e-15 1.325e-15
## Proportion of Variance 0.04212 0.01737 0.000e+00 0.000e+00 0.000e+00 0.000e+00
## Cumulative Proportion  0.98263 1.00000 1.000e+00 1.000e+00 1.000e+00 1.000e+00
##                             PC14
## Standard deviation     8.187e-16
## Proportion of Variance 0.000e+00
## Cumulative Proportion  1.000e+00

Observamos que con PC6 ya llegamos al 88% de la proporción de varianza explicada pero vemos algunos saltos interesantes. Para ver cuantas componentes explican el modelo hacemos la regla de codo, método heurístico para determinar el número óptimo de clusters en un conjunto de datos.

Al observar el gráfico vemos que la forma del codo sería en la 6|7 aunque el gráfico no es muy claro.

summary_pca <- summary(pca_result)
var_exp <- summary_pca$importance[2,]

barplot(var_exp, col = "#0BB363",
        main = "Gráfico PCA - Variance explained",
        xlab = "Principal Components",
        ylab = "Proportion of Variance Explained")

lines(var_exp, col = "#DC143C", type = "l", lwd = 2)

points(var_exp, col = "#DC143C", pch = 20, cex = 1.2)

A modo de ejercicio para comprender el funcionamiento de PCA puede resultar útil hacer estos cálculos, sin embargo para el dataset en cuestión ya se pueden formar modelos robustos con muy pocas variables (los m2 de superficie, por ejemplo, ya son suficientes para predecir un gran porcentaje de precios).

10 Aprendizaje supervisado

10.1 GLM

En aprendizaje automático, GLM (Generalized Linear Models) se refiere a una clase de modelos lineales que se utilizan para predecir una variable respuesta a partir de una o más variables explicativas. Tiene tres componentes: Variable respuesta, función de enlace y función lineal sistemática, y son muy útiles para utilizarse en clasificación, regresión y problemas de predicción.

Aplicado a nuestro dataset, lo utilizaremos para predecir el precio binario de las casas (cara/barata).

train_glm <- mhp_train[,c(3:7,10,11,12)]
test_glm <- mhp_test[,c(3:7,10,11,12)]
train_glm[,2:3]<- scale(train_glm[,2:3])
test_glm[,2:3]<- scale(test_glm[,2:3])

Tras seleccionar las variables y normalizarlas para que todas tengan el mismo peso, desarrollamos el modelo con la variable de respuesta “precio_bin” y las variables independientes “exterior”, “garage”, “elevator”, “rooms”, “precio_distrito” y “tipo_casa”, usando una distribución binomial.

glm_mhp = glm(formula = precio_bin ~ exterior + garage + elevator + rooms + m2 + precio_distrito + tipo_casa,
                 data = train_glm, 
                 family = binomial)
summary(glm_mhp)
## 
## Call:
## glm(formula = precio_bin ~ exterior + garage + elevator + rooms + 
##     m2 + precio_distrito + tipo_casa, family = binomial, data = train_glm)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.5396  -0.1705   0.0219   0.2836   5.0687  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.99261    0.62333  -1.592  0.11129    
## exterior0            -0.35269    0.13261  -2.660  0.00782 ** 
## garage0               0.79113    0.10737   7.368 1.73e-13 ***
## elevator0             1.67266    0.12813  13.054  < 2e-16 ***
## rooms                -0.01548    0.08151  -0.190  0.84940    
## m2                   -6.20689    0.23649 -26.246  < 2e-16 ***
## precio_distritomedio -2.52394    0.11266 -22.403  < 2e-16 ***
## precio_distritocaro  -4.33905    0.14476 -29.974  < 2e-16 ***
## tipo_casapiso         0.34751    0.62328   0.558  0.57716    
## tipo_casasotano       1.01944    0.62491   1.631  0.10282    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10718.8  on 7731  degrees of freedom
## Residual deviance:  3535.9  on 7722  degrees of freedom
## AIC: 3555.9
## 
## Number of Fisher Scoring iterations: 8

Tras realizar muchas comprobaciones y probar de maneras diversa el modelo no conseguimos interpretar bien los datos. No somos capaces de entender por qué sale inversamente proporcional los m2,el precio distrito mdio y el caro, que precisamente son las 3 variables que más definen que una casa sea cara, entonces no entendemos qué significa en este caso que salgan negativas. Lo dejamos pendiente para que por favor nos expliqueis a qué se puede deber esto y así saberlo para futuras ocasiones.

head(predict(glm_mhp))
##          1          2          3          4          5          6 
## -2.6264925  1.5662145 -6.1858748  0.1476553  0.4424557 -3.7636044

Esto sirve para predecir la probabilidad logit ajustado previamente y para mostrar las primeras seis predicciones. Cada número representa la probabilidad logit de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.

ajustados <- fitted(glm_mhp)
head(fitted(glm_mhp))
##           1           2           3           4           5           6 
## 0.067452746 0.827243283 0.002054072 0.536846908 0.608844012 0.022673932

Aquí obtenemos las probabilidades ajustadas que representan las probabilidades estimadas de que la variable dependiente “precio_bin” sea igual a 1 para cada observación en el conjunto de datos.

Gráficos de residuales: Estos gráficos permiten evaluar la calidad del ajuste del modelo. Puedes trazar residuales de Pearson, deviance, o residuales estandarizados frente a las predicciones ajustadas o valores ajustados para examinar si hay patrones en los residuales. Un buen ajuste del modelo no mostrará patrones claros en los residuales.

residuales <- residuals(glm_mhp, type = "pearson")
head(residuales)
##           1           2           3           4           5           6 
## -0.26894557 -2.18826117 -0.04536849 -1.07662124  0.80153404 -0.15231535
data_residuales <- data.frame(Residuales = residuales, Ajustados = ajustados)
ggplot(data_residuales, aes(x = Ajustados, y = Residuales)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(title = "Gráfico de residuales", x = "Valores ajustados", y = "Residuales") + ylim(-10,10)

predicciones_glm <- ifelse(test = glm_mhp$fitted.values > 0.5, no = 0, yes = 1)
matriz_confusion_glm <- table(glm_mhp$model$precio_bin, predicciones_glm,
                          dnn = c("observaciones", "predicciones"))
matriz_confusion_glm
##              predicciones
## observaciones    0    1
##             1 3449  414
##             0  329 3540

Finalmente generamos predicciones binarias basadas en las probabilidades ajustadas del modelo lineal generalizado y creamos una matriz de confusión para evaluar el rendimiento del modelo donde se muestra la distribución de las predicciones correctas e incorrectas en función de las observaciones reales y las predichas.

Obersvamos que el modelo no acierta nada, esto imagino que estará relacionado con lo anterior que no llegamos a entender.

La matriz se lee de la siguiente manera:

  • Verdaderos negativos (VN): La esquina superior izquierda representa el número de observaciones clasificadas correctamente como “0” (clase negativa).

  • Falsos positivos (FP): La esquina superior derecha representa el número de observaciones que en realidad eran “0” pero fueron clasificadas incorrectamente como “1” (clase positiva) por el modelo.

  • Falsos negativos (FN): La esquina inferior izquierda representa el número de observaciones que en realidad eran “1” pero fueron clasificadas incorrectamente como “0” por el modelo.

  • Verdaderos positivos (VP): La esquina inferior derecha representa el número de observaciones clasificadas correctamente como “1”.

10.2 KNN

El modelo KNN (K-Nearest Neighbors) es un algoritmo de aprendizaje supervisado que se utiliza tanto para clasificación como para regresión. Es un método basado en instancias que no hace supuestos sobre la forma funcional de la relación entre las variables independientes y la variable dependiente. KNN utiliza la información de los vecinos más cercanos para realizar predicciones. Para ello calcula distancias como la Euclidiana, Manhattan, Minkowski, etc…

La idea detrás de KNN es simple, una observación se clasifica o se le asigna un valor basado en las características de sus K vecinos más cercanos.

Algunas ventajas del algoritmo KNN son su simplicidad o su capacidad para manejar relaciones no lineales entre variables. Sin embargo también tiene algunas desventajas, como su sensibilidad a la maldición de la dimensionalidad, el efecto de las variables irrelevantes y el costo computacional asociado con el cálculo de distancias en conjuntos de datos grandes.

Procedemos al cálculo de este modelo e intentamos ver si hay un número correcto de k-vecinos que podemos usar.

long = 15
accuracy = rep(0,long)
f1score = rep(0,long)
recall = rep(0,long)
precision = rep(0,long)
for (i in 1:long)
{
  prediccion_knn_cv =knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], 
                            k=i, cl=mhp_train$precio_bin)
  accuracy[i] = sum(prediccion_knn_cv == mhp_train$precio_bin) /nrow(mhp_train)
  recall[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & mhp_train$precio_bin == TRUE) / sum(mhp_train$precio_bin == TRUE)
  precision[i] = sum(prediccion_knn_cv == mhp_train$precio_bin & prediccion_knn_cv == TRUE) / sum(prediccion_knn_cv == TRUE)
  f1score[i] = 2*precision[i]*recall[i]/(precision[i]+recall[i])
}
resultados_knn = as.data.frame(cbind(accuracy,f1score,precision,recall))
resultados_knn = resultados_knn %>% mutate(index=as.factor(seq(1:long)))

max(resultados_knn$f1score)
## [1] NaN
which.max(resultados_knn$f1score)
## integer(0)
ggplot(data=resultados_knn,aes(x=index,y=accuracy)) + 
  geom_col(colour="cyan4",fill="cyan3")+
  ggtitle("Accuracy")

ggplot(data=resultados_knn,aes(x=index,y=f1score)) + 
  geom_col(colour="orange4",fill="orange3") +
  ggtitle("F1_score values")

Tratamos de encontrar el valor óptimo de k utilizando validación cruzada. Se calculan varias métricas de evaluación, como precisión, recall y F1-score, para cada valor de k y finalmente tratamos de graficar la precisión en función del valor de k, pero analizando el gráfico no observamos ninguna diferencia, así que cogemos 5, que es el tamaño por convención que se suele coger.

Pasamos a realizar predicciones en el conjunto de datos de entrenamiento y en el de test además de realizar la matriz de confusión.

# En train
prediccion_knn5_train <- knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], 
                              k=5, cl=mhp_train$precio_bin, prob = TRUE)
confusionMatrix(table(prediccion_knn5_train,mhp_train$precio_bin), positive= "1")
## Confusion Matrix and Statistics
## 
##                      
## prediccion_knn5_train    1    0
##                     1 3233  531
##                     0  630 3338
##                                           
##                Accuracy : 0.8498          
##                  95% CI : (0.8417, 0.8577)
##     No Information Rate : 0.5004          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6997          
##                                           
##  Mcnemar's Test P-Value : 0.004026        
##                                           
##             Sensitivity : 0.8369          
##             Specificity : 0.8628          
##          Pos Pred Value : 0.8589          
##          Neg Pred Value : 0.8412          
##              Prevalence : 0.4996          
##          Detection Rate : 0.4181          
##    Detection Prevalence : 0.4868          
##       Balanced Accuracy : 0.8498          
##                                           
##        'Positive' Class : 1               
## 
#En test  
prediccion_knn5_test <- knn(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], mhp_test[,c("exterior","rooms","m2","elevator", "garage")],
                         k=5, cl=mhp_train$precio_bin, prob = TRUE)
confusionMatrix(table(prediccion_knn5_test,mhp_test$precio_bin), positive= "1")
## Confusion Matrix and Statistics
## 
##                     
## prediccion_knn5_test    1    0
##                    1 1577  250
##                    0  356 1683
##                                           
##                Accuracy : 0.8432          
##                  95% CI : (0.8314, 0.8546)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.6865          
##                                           
##  Mcnemar's Test P-Value : 1.996e-05       
##                                           
##             Sensitivity : 0.8158          
##             Specificity : 0.8707          
##          Pos Pred Value : 0.8632          
##          Neg Pred Value : 0.8254          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4079          
##    Detection Prevalence : 0.4726          
##       Balanced Accuracy : 0.8432          
##                                           
##        'Positive' Class : 1               
## 

La matriz de confusión en train nos muestra la cantidad de predicciones correctas e incorrectas de cada clase (1 y 0) en el conjunto de entrenamiento.

  • Accuracy (exactitud): 0.8482. Es la proporción de predicciones correctas en el conjunto de entrenamiento. En otras palabras, el modelo KNN clasifica correctamente el 84.82% de las observaciones del conjunto de entrenamiento.

  • Kappa: 0.6963. El índice kappa mide la concordancia entre las predicciones y los valores reales. Un kappa de 1 indica una concordancia perfecta, mientras que un kappa de 0 indica una concordancia totalmente aleatoria. Nuestro kappa de 0.6963 muestra una buena concordancia entre las predicciones y los valores reales en el conjunto de entrenamiento.

  • Mcnemar’s Test P-Value: 0.003201. Esta prueba evalúa si hay una diferencia significativa entre las predicciones falsas positivas y las predicciones falsas negativas.

La matriz de confusión en test nos muestra valores muy parecidos en los indicadores anteriores lo que nos lleva a pensar que el modelo tiene un rendimiento consistente y generaliza bien lo datos ya que el modelo clasifica correctamente las observaciones y no parece sobreajustado.

Hemos realizado diferentes pruebas quitando variables y se observa que si quitamos la variable m2 nos dice que hay muchos empates. Sin embargo, quitando las otras variables el acierto en train y test varía.

Obteniendo más o menos un 84% de acierto creemos que nos puede servir para clasificar las nuevas casas que entrasen en el dataset, puesto que no hemos sido capaces de mejorar ese %.

10.3 Decision Trees

El modelo de Árbol de Decisión es un algoritmo de aprendizaje supervisado que funciona dividiendo iterativamente el conjunto de datos en subconjuntos basados en las características del mismo. Al final de este proceso, se genera un árbol con nodos y hojas que representan las divisiones y las decisiones.

Algunas características importantes de los árboles de decisión son la interpretabilidad, capacidad para manejar datos no lineales, el poco preprocesamiento de los datos, la sensibilidad a datos desequilibrados y el posible sobreajuste en el caso de hacer un arbol con demasiadas ramificaciones.

A continuación hacemos el modelo de arbol de decisión.

arbol <- rpart(precio_bin ~ garage + elevator + rooms + tipo_casa + precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")

Para evaluar el modelo utilizamos la función predict sobre el conjunto test usando el argumento type = “class” que indica que se espera una salida categórica creandose una tabla de contingencias “tab1” con las predicciones que coinciden o no respecto a las observaciones reales.

Calculamos y obtenemos el número 0.8396275 que representa la tasa de acierto (accuracy) del modelo. El Árbol de Decisión hizo predicciones correctas en casi el 84% de los casos en el conjunto de datos de prueba. Nos parece un porcentaje bastante aceptable.

tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1
## [1] 0.8396275

Siguiendo con la validación del modelo podemos usar las funciones printcp y plotcp que se utilizan para analizar la complejidad y el rendimiento usando el “pruning” o poda de ramas del arbol que a menudo reduce su complejidad y permite una mejor generalización del mismo disminuyendo el subreajuste.

La función printcp es una función que muestra una tabla que contiene información sobre la complejidad y la tasa de error en diferentes niveles de tamaño del árbol proporcionando una visión general de cómo cambia el rendimiento a medida que aumenta la complejidad. El valor de cp ha de ser tal que la tasa de error de validación cruzada sea mínima.

La función plotcp crea un gráfico de la tasa de error. Este gráfico es útil para visualizar cómo se comporta el rendimiento del árbol a medida que aumenta su tamaño y ayuda a decidir qué tamaño de árbol es el más apropiado para el problema en cuestión.

printcp(arbol)
## 
## Classification tree:
## rpart(formula = precio_bin ~ garage + elevator + rooms + tipo_casa + 
##     precio_distrito, data = mhp_train, control = rpart.control(minsplit = 1))
## 
## Variables actually used in tree construction:
## [1] elevator        precio_distrito rooms          
## 
## Root node error: 3863/7732 = 0.49961
## 
## n= 7732 
## 
##         CP nsplit rel error  xerror      xstd
## 1 0.516438      0   1.00000 1.03055 0.0113762
## 2 0.044784      1   0.48356 0.48356 0.0097435
## 3 0.041677      3   0.39399 0.40098 0.0091108
## 4 0.016179      4   0.35232 0.35232 0.0086689
## 5 0.010000      6   0.31996 0.31996 0.0083418

Obtenidos los resultados, pasamos a analizarlos:

  • CP (Cost Complexity): Mide la complejidad asociado con cada nivel de poda en el árbol de decisión. Un CP más pequeño indica una complejidad menor.
  • nsplit: Número de divisiones (splits) realizadas hasta el momento en el árbol.
  • rel error: Error relativo en cada etapa de poda del árbol.
  • xerror: Error de validación cruzada (cross-validated) estimado para el árbol en cada etapa de poda.
  • xstd: Desviación estándar del error de validación cruzada (xerror).

En base a los resultados presentados, se puede observar que el árbol sin poda tiene un error relativo de 1.0 y un error de validación cruzada de 1.0189. A medida que se realizan más divisiones en el árbol el CP disminuye y también lo hace el error relativo y el error de validación cruzada.

Para seleccionar el árbol óptimo, generalmente se busca el CP más pequeño que tenga un error de validación cruzada dentro de una desviación estándar del mínimo error de validación cruzada.

arbol$cptable[which.min(arbol$cptable[, "xerror"]),
"CP"]
## [1] 0.01

El propósito es identificar el mejor árbol de decisión con respecto a su capacidad de generalización en datos no vistos. La función which.min() devuelve el índice del valor mínimo en la columna “xerror” de la tabla de complejidad cptable del árbol de decisión. Luego, con ese índice, se extrae el valor correspondiente de la columna “CP” en la tabla de complejidad.

plotcp(arbol)

Podamos el árbol como indica la función.

pruneTREE1 = prune(arbol, cp = arbol$cptable[which.min(arbol$cptable[,
"xerror"]), "CP"])
fancyRpartPlot(pruneTREE1, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")

pruneTREE2 = prune(arbol, cp = arbol$cptable[5, "CP"])
fancyRpartPlot(pruneTREE2, uniform = TRUE, main = "Pruned Classification Tree",
sub = "")

tab2 = table(pred = predict(pruneTREE2, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
acierto2 = sum(diag(tab2))/ntest
acierto2
## [1] 0.8396275

Tras la poda nos sale literalmente es mismo resultado.

Creamos la matriz de confusión.

prediccion_1 <- predict(arbol, newdata = mhp_train, type = "class")
confusionMatrix(prediccion_1, mhp_train[["precio_bin"]])
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    0
##          1 3246  619
##          0  617 3250
##                                           
##                Accuracy : 0.8401          
##                  95% CI : (0.8318, 0.8482)
##     No Information Rate : 0.5004          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.6803          
##                                           
##  Mcnemar's Test P-Value : 0.9773          
##                                           
##             Sensitivity : 0.8403          
##             Specificity : 0.8400          
##          Pos Pred Value : 0.8398          
##          Neg Pred Value : 0.8404          
##              Prevalence : 0.4996          
##          Detection Rate : 0.4198          
##    Detection Prevalence : 0.4999          
##       Balanced Accuracy : 0.8401          
##                                           
##        'Positive' Class : 1               
## 

Los resultados en la matriz de confusión son los siguientes:

  • Verdaderos positivos (VP): 3246
  • Verdaderos negativos (VN): 3250
  • Falsos positivos (FP): 619
  • Falsos negativos (FN): 617

También obtenemos:

  • Accuracy (exactitud): 84.01% de predicciones correctas.
  • Kappa: 68.03%
  • Sensitivity (sensibilidad): 0.8403. También conocida como recall o tasa verdadera positiva, mide la proporción de casos positivos que se identificaron correctamente.
  • Specificity (especificidad): 0.8400. Mide la proporción de casos negativos que se identificaron correctamente. Un valor más alto es mejor.
  • Pos Pred Value (valor predictivo positivo): 0.8398. Mide la proporción de predicciones positivas que son realmente positivas. Un valor más alto es mejor.
  • Neg Pred Value (valor predictivo negativo): 0.8404. Mide la proporción de predicciones negativas que son realmente negativas. Un valor más alto es mejor.
  • Mcnemar’s Test P-Value: 0.9773. Esta prueba compara las diferencias entre los falsos positivos y los falsos negativos. Un valor alto de “p” (como en este caso) sugiere que no hay una diferencia significativa entre ellos.

En resumen, este modelo de árbol de decisión podado tiene una exactitud del 84.01% y un coeficiente kappa de 0.6803, lo que indica un buen rendimiento en la clasificación.

Sin embargo para este modelo se han tenido en cuenta bastantes variables, y aunque el rendimiento es bueno, se puede alcanzar un rendimiento casi igual utilizando sólo los metros cuadrados, como se puede ver a continuación

arbol <- rpart(precio_bin ~ m2, data = mhp_train, control = rpart.control(minsplit = 1))
fancyRpartPlot(arbol, sub = "")

tab1 = table(pred = predict(arbol, mhp_test, type = "class"),
obs = mhp_test$precio_bin)
ntest = nrow(mhp_test)
acierto1 = sum(diag(tab1))/ntest
acierto1
## [1] 0.8344542

Esta es una forma simple de confirmar la relevancia que tiene la superficie de una casa sobre su precio final, en contraste con otras variables como el tener ascensor o garaje. Además, esto nos deja un modelo mucho más simple y explicable.

10.4 Random Forest

Algoritmo de aprendizaje supervisado que se basa en la construcción de múltiples árboles de decisión durante el entrenamiento y la combinación de sus resultados para hacer una predicción. El objetivo principal de este enfoque es mejorar la precisión y la robustez del modelo en comparación con un único árbol de decisión.

El funcionamiento de un algoritmo Random Forest sería el siguiente:

  • Seleccionar muestras de arranque: El algoritmo selecciona aleatoriamente un subconjunto de muestras del conjunto de datos de entrenamiento, con reemplazo. Esto significa que una muestra puede ser seleccionada más de una vez.

  • Crear un árbol de decisión: Se crea un árbol de decisión utilizando el subconjunto de muestras seleccionado. Durante la construcción del árbol, en cada división del nodo, se selecciona un subconjunto aleatorio de características como candidatas para la división. La mejor división se elige en función de la medida de impureza como la entropía o el índice de Gini. Esto introduce aleatoriedad adicional en el proceso, lo que ayuda a reducir la correlación entre los árboles.

  • Repetir: El proceso se repite varias veces (número de árboles especificado) para construir un conjunto de árboles de decisión.

  • Combinar las predicciones: Para hacer una predicción en un nuevo caso, el algoritmo Random Forest ejecuta todos los árboles individualmente y luego combina sus resultados. En el caso de la clasificación, la clase con más votos (predicciones) de todos los árboles se elige como la predicción final.

Pasamos a construir el modelo comenzando con 10 árboles aleatórios.

# Escalamos la columna 4 en train y test
set.seed(19042023)
mhp_train_scaled <- mhp_train
mhp_train_scaled[, 4] <- scale(mhp_train[, 4])

mhp_test_scaled <- mhp_test
mhp_test_scaled[, 4] <- scale(mhp_test[, 4])

# Nos quedamos con todas las variables menos barrio y distrito porque hay muchas categorias
# Cogemos solo 10 árboles en el parámetro ntree
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
                           y = mhp_train_scaled$precio_bin,
                           ntree = 10)
# Predicción de los resultados con el conjunto de testing
y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])
confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    0
##          1 1676  127
##          0  257 1806
##                                           
##                Accuracy : 0.9007          
##                  95% CI : (0.8908, 0.9099)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8013          
##                                           
##  Mcnemar's Test P-Value : 4.61e-11        
##                                           
##             Sensitivity : 0.8670          
##             Specificity : 0.9343          
##          Pos Pred Value : 0.9296          
##          Neg Pred Value : 0.8754          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4335          
##    Detection Prevalence : 0.4664          
##       Balanced Accuracy : 0.9007          
##                                           
##        'Positive' Class : 1               
## 

Ahora probamos con 100 árboles aleatorios para ver si cambia mucho el resultado.

set.seed(19042022)
classifier <- randomForest(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
                           y = mhp_train_scaled$precio_bin,
                           ntree = 100)
y_pred = predict(classifier, newdata = mhp_test_scaled[,c(3,4,5,6,7,10,11)])
confusionMatrix(y_pred, mhp_test_scaled[["precio_bin"]])
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    1    0
##          1 1684  108
##          0  249 1825
##                                           
##                Accuracy : 0.9077          
##                  95% CI : (0.8981, 0.9166)
##     No Information Rate : 0.5             
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.8153          
##                                           
##  Mcnemar's Test P-Value : 1.267e-13       
##                                           
##             Sensitivity : 0.8712          
##             Specificity : 0.9441          
##          Pos Pred Value : 0.9397          
##          Neg Pred Value : 0.8799          
##              Prevalence : 0.5000          
##          Detection Rate : 0.4356          
##    Detection Prevalence : 0.4635          
##       Balanced Accuracy : 0.9077          
##                                           
##        'Positive' Class : 1               
## 

La precisión del modelo es del 0,9077 y el intervalo de confianza del 95% es (0,8908, 0,9099), lo que indica que podemos estar bastante seguros de que la precisión real del modelo está dentro de este rango. El valor de kappa de 0,8013 indica que hay una buena concordancia entre las predicciones del modelo y las observaciones reales.

Además, la sensibilidad del modelo es del 0,8670 y la especificidad del 0,9343. La tasa de detección es del 0,4335 y la prevalencia es del 0,5000. En general, el modelo tiene un desempeño bastante bueno en la clasificación de las viviendas por encima y por debajo del umbral.

10.5 SVM

Support Vector Machine, algoritmo de aprendizaje supervisado que se utiliza tanto para problemas de clasificación como de regresión. En el caso de la clasificación, el objetivo principal de SVM es encontrar el hiperplano óptimo que separa los datos en diferentes clases. Un hiperplano es un límite de decisión que divide el espacio de características en dos partes, de modo que cada parte contenga puntos de datos de una clase específica. En el caso de la regresión, SVM busca encontrar una función que tenga el menor error de ajuste posible.

El algoritmo SVM funciona de la siguiente manera:

  • Encuentra el hiperplano óptimo: El algoritmo busca el hiperplano que tenga el margen máximo entre las dos clases. El margen es la distancia entre el hiperplano y los puntos de datos más cercanos a él, llamados “vectores de soporte”. Estos vectores de soporte son los puntos de datos que definen el límite de decisión y son fundamentales para el proceso de SVM.

  • Solución de un problema de optimización: La búsqueda del hiperplano óptimo implica la solución de un problema de optimización convexa. Este problema se resuelve mediante técnicas de optimización como el método de Lagrange o el algoritmo de optimización secuencial mínima (SMO).

  • Transformación a través del kernel: En casos en los que los datos no son linealmente separables en el espacio de características original, SVM puede utilizar una función de kernel para transformar los datos en un espacio de características de mayor dimensión donde se puedan separar linealmente. Los kernels comunes incluyen el kernel lineal, el kernel polinómico, el kernel de base radial (RBF) y el kernel sigmoide. La elección del kernel adecuado es crucial para el rendimiento del modelo SVM.

  • Clasificación: Para hacer una predicción en un nuevo punto de datos, el algoritmo SVM evalúa la posición del punto con respecto al hiperplano óptimo. Si el punto está por encima del hiperplano, se asigna a una clase, y si está por debajo, se asigna a la otra clase.

Para resolver SVM se necesita conocer el kernel y el parámetro de regularización, que es una cuota superior sobre las variables duales, en definitiva, un coste. Pasamos a construirun primer el modelo usando un kernel radial RBF.

# Estandarizamos las variables numéricas "m2" y "rooms"
mhp_train_svm <- mhp_train
mhp_train_svm$m2 <- scale(mhp_train$m2)
mhp_train_svm$rooms <- scale(mhp_train$rooms)

mhp_test_svm <- mhp_test
mhp_test_svm$m2 <- scale(mhp_test$m2)
mhp_test_svm$rooms <- scale(mhp_test$rooms)

# Entrenamos el modelo SVM  kernel con kernel radial

modelo_svm <- svm(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito, data = mhp_train_svm, kernel = "radial", cost = 10, scale = TRUE)

# Hacemos predicciones en el conjunto de prueba
svm_predicciones <- predict(modelo_svm, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")])


# Calculamos la matriz de confusión y la precisión del modelo
svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test_svm$precio_bin)
print(svm_matriz_confusion)
##           Real
## Prediccion    1    0
##          1 1692  116
##          0  241 1817
precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
print(paste0("Precisión del modelo SVM: ", round(precision * 100, 2), "%"))
## [1] "Precisión del modelo SVM: 90.77%"

Esto significa que:

  • 1692 casas con precio alto fueron clasificadas correctamente como precio alto (verdaderos positivos).
  • 1817 casas con precio bajo fueron clasificadas correctamente como precio bajo (verdaderos negativos).
  • 241 casas con precio alto fueron clasificadas incorrectamente como precio bajo (falsos negativos).
  • 116 casas con precio bajo fueron clasificadas incorrectamente como precio alto (falsos positivos).

La precisión del modelo es del 90.77%.

Elección del Kernel

La elección del kernel en un modelo SVM depende de la naturaleza y la distribución de los datos. No hay una regla única para seleccionar el kernel, pero aquí hay algunas pautas generales para ayudarte a tomar una decisión:

  • Kernel lineal: Si tus datos son linealmente separables o la cantidad de atributos es alta en comparación con la cantidad de muestras, entonces un kernel lineal suele funcionar bien. Además, es computacionalmente más eficiente que otros kernels.

  • Kernel polinomial: Si tus datos tienen una estructura que puede ser capturada por un polinomio de un grado específico, el kernel polinomial podría ser una buena opción. Sin embargo, a medida que aumenta el grado del polinomio, el riesgo de sobreajuste aumenta y la complejidad computacional puede crecer.

  • Kernel radial (RBF): El kernel radial es una opción popular y versátil en muchos casos, especialmente cuando no se sabe de antemano si los datos son linealmente separables. Puede manejar datos no lineales y complejos. Sin embargo, encontrar los parámetros adecuados (costo y gamma) puede requerir una búsqueda en el espacio de parámetros y, en consecuencia, un mayor tiempo de cómputo.

  • Kernel sigmoide: El kernel sigmoide es similar a una función de activación de una red neuronal y puede ser una opción si deseas un enfoque similar al de una red neuronal pero con un modelo SVM.

Vamos a probar a representar los datos para ver si de alguna manera puede tomarse una mejor decisión sobre el Kernel a usar.

mhp_train %>%
  ggplot(aes(x = m2, y = rooms, color = as.factor(precio_bin))) +
  geom_point() +
  facet_grid(elevator ~ garage, labeller = labeller(
    elevator = c('0' = 'Sin Elevador', '1' = 'Con Elevador'),
    garage = c('0' = 'Sin Garaje', '1' = 'Con Garaje')
  )) +
  labs(title = "Diagrama de dispersión de m2 vs rooms, según garage y elevator",
       x = "Metros cuadrados",
       y = "Número de habitaciones",
       color = "Precio bin") +
  scale_color_manual(values = c("1" = "#0BB363", "0" = "#DC143C")) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"),
        plot.title.position = "plot")

Probamos varios Kernels.

Creamos una función entrenar_evaluar_svm() que entrenará y evaluará un modelo SVM utilizando el tipo de kernel especificado en el argumento kernel_type. Luego con lapply() se aplicará esta función a cada tipo de kernel en el vector kernels. La función presentará un matriz de confusión y la precisión para cada tipo de kernel.

# Definir una función para entrenar y evaluar un modelo SVM con un tipo de kernel específico
entrenar_evaluar_svm <- function(kernel_type) {
  # Entrenamos el modelo SVM con el kernel especificado
  modelo_svm <- svm(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito, data = mhp_train_svm, kernel = kernel_type, cost = 10, scale = TRUE)

  # Hacemos predicciones en el conjunto de prueba
  svm_predicciones <- predict(modelo_svm, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")])

  # Calculamos la matriz de confusión y la precisión del modelo
  svm_matriz_confusion <- table(Prediccion = svm_predicciones, Real = mhp_test_svm$precio_bin)
  precision <- sum(diag(svm_matriz_confusion)) / sum(svm_matriz_confusion)
  
  # Devolvemos una lista con la precisión y la matriz de confusión
  return(list(precision = precision, matriz_confusion = svm_matriz_confusion))
}

# Probar diferentes tipos de kernel
kernels <- c("linear", "polynomial", "radial", "sigmoid")
resultados <- lapply(kernels, entrenar_evaluar_svm)

# Imprimir los resultados
for (i in 1:length(kernels)) {
  cat(paste("Precisión del modelo SVM con kernel", kernels[i], ":", round(resultados[[i]]$precision * 100, 2), "%\n"))
  print(resultados[[i]]$matriz_confusion)
}
## Precisión del modelo SVM con kernel linear : 90.46 %
##           Real
## Prediccion    1    0
##          1 1706  142
##          0  227 1791
## Precisión del modelo SVM con kernel polynomial : 90.43 %
##           Real
## Prediccion    1    0
##          1 1688  125
##          0  245 1808
## Precisión del modelo SVM con kernel radial : 90.77 %
##           Real
## Prediccion    1    0
##          1 1692  116
##          0  241 1817
## Precisión del modelo SVM con kernel sigmoid : 85 %
##           Real
## Prediccion    1    0
##          1 1633  280
##          0  300 1653

Se observa que el KERNEL RADIAL es el que mejor predice, ahora vamos a ajustar algunos parámetros del modelo utilizando la función tune para encontrar los mejores parámetros.

# Entrenar un modelo SVM con diferentes parámetros
modelo_svm_radial <- svm(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito, data = mhp_train_svm, kernel = "radial", cost = 10, scale = TRUE)

# Utilizar la función tune() para encontrar los mejores parámetros
tuned_parameters <- tune(svm, 
                         precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito,
                         data = mhp_train_svm,
                         kernel = "radial",
                         ranges = list(cost = c(0.1, 1, 10, 100, 1000),
                                       gamma = c(0.1, 0.5, 1, 2, 5)))

# Entrenar un modelo SVM con los mejores parámetros encontrados
best_model <- tuned_parameters$best.model

# Hacer predicciones en el conjunto de prueba con el mejor modelo
best_predictions <- predict(best_model, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")])

# Calcular la matriz de confusión y la precisión del mejor modelo
best_matriz_confusion <- table(Prediccion = best_predictions, Real = mhp_test_svm$precio_bin)
print(best_matriz_confusion)
##           Real
## Prediccion    1    0
##          1 1692  118
##          0  241 1815
best_precision <- sum(diag(best_matriz_confusion)) / sum(best_matriz_confusion)
print(paste0("Precisión del mejor modelo SVM: ", round(best_precision * 100, 2), "%"))
## [1] "Precisión del mejor modelo SVM: 90.71%"

La precisión del modelo no mejora del 90.77% anterior y el coste computacional en la busqueda del mejor modelo es enorme. Parece que los parámetros del modelo elegidos durante el proceso de afinación no han sido los correctos. La afinación e hiperparámetros implica buscar la mejor combinación de hiperparámetros dentro del rango de búsqueda que se proporciona. Si el rango de búsqueda no contiene la combinación óptima de hiperparámetros, el proceso de afinación podría terminar eligiendo una combinación subóptima que resulte en un rendimiento inferior.

11 Evaluación y comparación de modelos

Evaluamos y comparamos los modelos que mayor porcentaje de predicción nos han proporcionado, Random Forest y SVM radial cost 10 (modelo inicial que nos daba un alto porcentaje de predicción y bajo coste computacional).

11.1 Random Forest

11.1.1 K-fold Cross Validation

K-fold Cross Validation es un método de validación cruzada donde la idea principal es dividir el conjunto de datos en K subconjuntos más pequeños, llamados “folds”, y utilizarlos para entrenar y probar el modelo de manera iterativa. Esto ayuda a evitar el sobreajuste y proporciona una estimación más precisa del rendimiento del modelo en datos no vistos.

El proceso de K-fold Cross Validation es el siguiente:

  • Dividir el conjunto de datos en K subconjuntos (folds) de aproximadamente el mismo tamaño.
    • Para cada fold:
        1. Utilizar K-1 folds para entrenar el modelo.
        1. Utilizar el fold restante como conjunto de prueba para evaluar el rendimiento del modelo.
  • Repetir el proceso K veces, asegurándose de que cada fold se utilice como conjunto de prueba exactamente una vez.
  • Calcular métricas de rendimiento, como precisión, error cuadrático medio, etc., para cada una de las K iteraciones y promediar los resultados para obtener una estimación más sólida del rendimiento del modelo.

La ventaja principal de utilizar K-fold Cross Validation es que todos los datos se utilizan tanto para entrenamiento como para pruebas, lo que proporciona una mejor estimación del rendimiento del modelo. Además, al repetir el proceso K veces, se reduce la varianza en las métricas de rendimiento, lo que resulta en una evaluación más confiable del modelo.

# Renombramos los niveles de la variable de respuesta en los conjuntos de datos de entrenamiento y prueba
mhp_train_scaled$precio_bin <- as.factor(mhp_train_scaled$precio_bin)
levels(mhp_train_scaled$precio_bin) <- c("Level0", "Level1")

mhp_test_scaled$precio_bin <- as.factor(mhp_test_scaled$precio_bin)
levels(mhp_test_scaled$precio_bin) <- c("Level0", "Level1")

# Configuramos la validación cruzada K-fold (por ejemplo, K=10)
k_folds <- 10
folds <- createFolds(mhp_train_scaled$precio_bin, k = k_folds)

# Ajustamos un modelo Random Forest usando validación cruzada K-fold con la librería caret
control <- trainControl(method = "cv", number = k_folds, index = folds)
K_F_modelo_rf <- train(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10, 11)],
                   y = mhp_train_scaled$precio_bin,
                   method = "rf",
                   ntree = 10,
                   trControl = control)

# Imprimimos el resultado del modelo
print(K_F_modelo_rf)
## Random Forest 
## 
## 7732 samples
##    7 predictor
##    2 classes: 'Level0', 'Level1' 
## 
## No pre-processing
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 773, 773, 772, 774, 773, 773, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.8911737  0.7823444
##   4     0.8847361  0.7694716
##   7     0.8757834  0.7515664
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

Estos resultados muestran el rendimiento del modelo de Random Forest para diferentes valores del hiperparámetro mtry que es el número de predictores que se utilizan en cada árbol de decisión. Se han probado 3 valores diferentes: 2, 4 y 7. y para cada valor de mtry se muestran precisiones cercanas al 90% como en el modelo inicial. K-fold validation nos indica que el mayor accuracy se obtiene con tan solo dos predictores, no hay necesidad de más.

11.1.2 N-fold Cross Validation

Es otra técnica de validación donde el conjunto de datos se divide en N subconjuntos y se entrena N veces utilizando cada vez un pliegue diferente como conjunto de validación, y los restantes N-1 pliegues como conjunto de entrenamiento.

Al final del proceso, se obtienen N métricas de rendimiento, una por cada iteración. Estas métricas se promedian para obtener una estimación única del rendimiento del modelo. Esta estimación es más fiable que la que se obtiene utilizando únicamente un único conjunto de entrenamiento y validación, ya que reduce el riesgo de que el rendimiento del modelo se vea afectado por una división aleatoria específica del conjunto de datos.

La validación cruzada de N pliegues ayuda a identificar si un modelo está sobreajustado (overfitting) o infraajustado (underfitting).

# Configura el método de control de entrenamiento
control <- trainControl(method = "cv",
                        number = 5, # N-fold Cross Validation
                        savePredictions = "final",
                        classProbs = TRUE)

# Entrena el modelo de Random Forest con N-fold Cross Validation
N_F_modelo_rf <- train(x = mhp_train_scaled[, c(3, 4, 5, 6, 7, 10,11)],
                   y = mhp_train_scaled$precio_bin,
                   method = "rf",
                   trControl = control,
                   tuneGrid = data.frame(mtry = c(2, 4, 7)))

# Muestra los resultados del modelo
print(N_F_modelo_rf)
## Random Forest 
## 
## 7732 samples
##    7 predictor
##    2 classes: 'Level0', 'Level1' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6186, 6185, 6186, 6186, 6185 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.9053296  0.8106515
##   4     0.9046822  0.8093612
##   7     0.8909726  0.7819452
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

Para hacer N-fold Validation se entrenó un modelo de Random Forest utilizando validación cruzada de 5 pliegues, y se encontró que el valor óptimo de mtry para este conjunto de datos y clasificación es 2. El modelo alcanzó una exactitud promedio del 90.4% y un índice Kappa promedio de 0.809.

11.2 SVM

Repetimos las validaciones para SVM.

11.2.1 K-fold Cross Validation

# fórmula
formula_svm <- precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito

# Configura los parámetros de entrenamiento
control_svm <- trainControl(method = "cv", number = 10) # 10-fold Cross Validation

# Entrena el modelo SVM con K-fold Cross Validation
K_F_modelo_svm_cv <- train(formula_svm,
                       data = mhp_train_svm,
                       method = "svmRadial",
                       trControl = control_svm,
                       preProcess = c("center", "scale"),
                       tuneGrid = data.frame(sigma = 0.1, C = 10))

# Muestra los resultados
print(K_F_modelo_svm_cv)
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 7732 samples
##    6 predictor
##    2 classes: '1', '0' 
## 
## Pre-processing: centered (8), scaled (8) 
## Resampling: Cross-Validated (10 fold) 
## Summary of sample sizes: 6960, 6959, 6958, 6959, 6959, 6959, ... 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.9027401  0.8054708
## 
## Tuning parameter 'sigma' was held constant at a value of 0.1
## Tuning
##  parameter 'C' was held constant at a value of 10

Le realizó la validación cruzada con 6 predictores y 10 pliegues alcanzándose una precisión del 90.43% y un Kappa del 80.86% mostrando una alta concordancia entre las predicciones y los valores reales ajustados por el azar.

11.2.2 N-fold Cross Validation

# Cambia los niveles de la variable de clase 'precio_bin' a nombres válidos
mhp_train_svm$precio_bin <- factor(mhp_train_svm$precio_bin, labels = c("Level0", "Level1"))

# Definimos los controles de entrenamiento para la validación cruzada N-fold
control <- trainControl(method = "cv", number = 5, savePredictions = "final", classProbs = TRUE)

# Creamos el modelo SVM con validación cruzada N-fold
N_F_modelo_svm_cv <- train(precio_bin ~ rooms + m2 + garage + elevator + tipo_casa + precio_distrito,
                       data = mhp_train_svm,
                       method = "svmRadial",
                       trControl = control,
                       preProcess = c("center", "scale"),
                       tuneGrid = data.frame(sigma = 0.1, C = 10))

# Mostramos los resultados del modelo SVM con validación cruzada N-fold
print(N_F_modelo_svm_cv)
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 7732 samples
##    6 predictor
##    2 classes: 'Level0', 'Level1' 
## 
## Pre-processing: centered (8), scaled (8) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6185, 6186, 6187, 6185, 6185 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.9032597  0.8065112
## 
## Tuning parameter 'sigma' was held constant at a value of 0.1
## Tuning
##  parameter 'C' was held constant at a value of 10
knn_model <- knn.cv(mhp_train[,c("exterior","rooms","m2","elevator", "garage")], 
                    k = 5, cl = mhp_train$precio_bin)

Le realizó la validación cruzada con 6 predictores y 10 pliegues alcanzándose una precisión del 90.29% y un Kappa del 80.57%.

12 Elección del mejor punto de corte según la curva ROC

Nos disponemos a calcular las curvas ROC para cada uno de los modelo y los visualizamos en un mismo plot para mejor comparabilidad.

# GLM
probs_glm <- predict(glm_mhp, test_glm, type = "response")

# KNN
probs_knn <- attr(prediccion_knn5_test, "prob")

# Decisión Tree
probs_dt <- predict(arbol, mhp_test, type = "prob")[, 2]

# Random Forest
probs_rf <- predict(classifier, mhp_test, type = "prob")[, 2]

# SVM
probs_svm <- predict(modelo_svm, mhp_test_svm[, c("rooms", "m2", "garage", "elevator", "tipo_casa", "precio_distrito")], probability = TRUE)
probs_svm <- attr(probs_svm, "probabilities")[, 2]
# Calcula las curvas ROC
roc_glm <- roc(mhp_test$precio_bin, probs_glm)
roc_knn <- roc(mhp_test$precio_bin, probs_knn)
roc_dt <- roc(mhp_test$precio_bin, probs_dt)
roc_rf <- roc(mhp_test$precio_bin, probs_rf)
#roc_svm <- roc(mhp_test$precio_bin, probs_svm)

# Configura el gráfico
plot(roc_glm, col = "blue", lwd = 2, legacy.axes = TRUE, main = "Curvas ROC")
lines(roc_knn, col = "red", lwd = 2)
lines(roc_dt, col = "green", lwd = 2)
lines(roc_rf, col = "darkred", lwd = 2)
#lines(roc_svm, col = "orange", lwd = 2)

# Añade leyenda
legend("bottomright", legend = c("GLM", "KNN", "Decisión Tree", "Random Forest", "SVM"),
       col = c("blue", "red", "green", "darkred", "orange"), lwd = 2)

CONCLUSIONES

Después de evaluar todos los modelos de aprendizaje automático, se observa que los mejores resultados se obtienen con el modelo de Random Forest y el modelo SVM (Support Vector Machine). Sin embargo, el modelo Random Forest es preferido sobre el modelo SVM por las siguientes razones:

  1. Coste computacional: El modelo Random Forest generalmente tiene un menor costo computacional en comparación con el modelo SVM, especialmente cuando se trabaja con grandes conjuntos de datos. Esto significa que el tiempo de entrenamiento y predicción es más corto y requiere menos recursos computacionales.

  2. Simplicidad del modelo: Aunque ambos modelos son eficaces, el modelo Random Forest es más fácil de entender y de explicar. Los árboles de decisión que forman el Random Forest pueden visualizarse y analizarse fácilmente, lo que facilita la interpretación de los resultados y la identificación de características importantes.

  3. Explicabilidad: La explicabilidad es un aspecto crítico en la ciencia de datos y el aprendizaje automático. Los modelos que son más fáciles de explicar y entender son preferibles en muchas situaciones, especialmente cuando se necesita justificar las decisiones basadas en estos modelos. El modelo Random Forest ofrece una mejor explicabilidad en comparación con el modelo SVM, ya que es más fácil de visualizar y entender cómo se toman las decisiones en el modelo.

  4. Rendimiento similar: A pesar de las ventajas mencionadas anteriormente, la elección entre Random Forest y SVM se basa en gran medida en su rendimiento en la tarea específica. En este caso, ambos modelos proporcionan resultados similares, lo que respalda la decisión de elegir el modelo Random Forest debido a sus otras ventajas.

13 Deep Learning

El Deep Learning es una rama del aprendizaje automático (machine learning) que se enfoca en el entrenamiento de redes neuronales artificiales profundas para aprender y extraer patrones complejos a partir de datos. Estas redes neuronales se componen de múltiples capas de nodos interconectados, lo que les permite aprender representaciones de datos en diferentes niveles de abstracción.

A diferencia de otros métodos de aprendizaje automático, el Deep Learning se caracteriza por su capacidad para aprender automáticamente características relevantes directamente de los datos, sin necesidad de que sean especificadas o extraídas manualmente. Esto lo logra mediante el proceso de entrenamiento, donde la red neuronal ajusta los pesos y los sesgos de sus conexiones para minimizar una función de pérdida y mejorar la capacidad de predicción o clasificación.

El Deep Learning ha mostrado un gran éxito en una amplia variedad de tareas de aprendizaje automático, como reconocimiento de imágenes, procesamiento de lenguaje natural, reconocimiento de voz, traducción automática… Esto se debe a su capacidad para modelar relaciones y patrones complejos en los datos, así como a la disponibilidad de grandes conjuntos de datos y avances en hardware de computación como los GPUs que aceleran el entrenamiento de redes neuronales profundas.

Algunas de las arquitecturas de redes neuronales profundas más utilizadas en Deep Learning incluyen las redes neuronales convolucionales (CNN) para el procesamiento de imágenes, redes neuronales recurrentes (RNN) para el procesamiento de secuencias y las redes neuronales generativas adversariales (GAN) para la generación de contenido nuevo.

Entonces… ¿Cómo funcionan?

Se podría decir que funcionan como un cerebro humano, hay neuronas y cada una de ellas recibe diferentes señales de entrada (𝑥1, 𝑥2, …) de fuentes externas o de otras neuronas. Cada señal de entrada se multiplica por un valor denominado peso que da una idea de la fuerza de dicha conexión (𝑤1, 𝑤2, …). La fase de entrenamiento modifica los pesos, y la neurona calcula una salida, usando Σ y una función de activación 𝑓.

Aplicado a nuestro dataset del mercado inmobiliario en Madrid, vamos a ir entrenando una red neuronal con distintas configuraciones de capas, neuronas y funciones de activación para ver si obtenemos un modelo realmente predictivo del precio de los inmuebles y por consiguiente generar una buena herramienta de clasificación de los mismos.

14 Explicabilidad

15 Series Temporales



Continuando con nuestro proyecto de Machine Learning pero apartándonos bastante del dataset original basado en la predicción del precio y etiquetado de casas en Madrid, procedemos a hacer un estudio de ventas en un famoso centro comercial estadounidense, Walmart.

Para ello hemos obtenido un dataset bastante completo con información de ventas de dicho establecimiento de un repositorio en la siguiente web (kaggle - Walmart Sales Dataset of 45 Stores) con formato long format.

Para el análisis de series temporales conviene saber que son un conjuntos de datos secuenciales que representan observaciones realizadas a lo largo del tiempo, generalmente en intervalos regulares. En este tipo de datos, el orden temporal es fundamental, ya que cada observación está asociada a un momento específico en el tiempo.

El análisis de series temporales se enfoca en comprender y modelar el patrón temporal subyacente en los datos, con el objetivo de predecir o explicar su comportamiento futuro. Las series temporales pueden presentar diferentes características, como tendencias, estacionalidad, ciclos y componentes aleatorios.

Al analizar una serie temporal, se busca identificar patrones y estructuras temporales, como cambios de tendencia, patrones estacionales o efectos de eventos específicos. Esto se logra mediante técnicas estadísticas y modelos matemáticos que permiten modelar y predecir el comportamiento futuro de la serie.

Las series temporales pueden utilizarse en una amplia gama de aplicaciones, como pronósticos económicos, análisis de ventas, previsión del clima, análisis de datos financieros…

En el contexto de nuestro proyecto con el dataset de ventas de Waltmark, nos disponemos a utilizar técnicas de análisis de series temporales para estudiar el comportamiento de las ventas a lo largo del tiempo, identificar patrones estacionales, detectar cambios en las tendencias, analizar el impacto de eventos especiales y realizar pronósticos de las ventas futuras.

15.1 Lectura y preparación de los datos

walmart <- read_csv("walmart.csv")
head(walmart, 10) %>%
  kbl() %>%
  kable_material(c("striped", "hover")) %>%
  scroll_box(width = "100%", height = "350px")
Store Date Weekly_Sales Holiday_Flag Temperature Fuel_Price CPI Unemployment
1 05-02-2010 1643691 0 42.31 2.572 211.0964 8.106
1 12-02-2010 1641957 1 38.51 2.548 211.2422 8.106
1 19-02-2010 1611968 0 39.93 2.514 211.2891 8.106
1 26-02-2010 1409728 0 46.63 2.561 211.3196 8.106
1 05-03-2010 1554807 0 46.50 2.625 211.3501 8.106
1 12-03-2010 1439542 0 57.79 2.667 211.3806 8.106
1 19-03-2010 1472516 0 54.58 2.720 211.2156 8.106
1 26-03-2010 1404430 0 51.45 2.732 211.0180 8.106
1 02-04-2010 1594968 0 62.27 2.719 210.8204 7.808
1 09-04-2010 1545419 0 65.86 2.770 210.6229 7.808
walmart$Store <- factor(walmart$Store)
walmart$Date <- as.Date(walmart$Date, format = "%d-%m-%Y")
walmart$Holiday_Flag <- as.factor(walmart$Holiday_Flag)
summary(walmart)
##      Store           Date             Weekly_Sales     Holiday_Flag
##  1      : 143   Min.   :2010-02-05   Min.   : 209986   0:5985      
##  2      : 143   1st Qu.:2010-10-08   1st Qu.: 553350   1: 450      
##  3      : 143   Median :2011-06-17   Median : 960746               
##  4      : 143   Mean   :2011-06-17   Mean   :1046965               
##  5      : 143   3rd Qu.:2012-02-24   3rd Qu.:1420159               
##  6      : 143   Max.   :2012-10-26   Max.   :3818686               
##  (Other):5577                                                      
##   Temperature       Fuel_Price         CPI         Unemployment   
##  Min.   : -2.06   Min.   :2.472   Min.   :126.1   Min.   : 3.879  
##  1st Qu.: 47.46   1st Qu.:2.933   1st Qu.:131.7   1st Qu.: 6.891  
##  Median : 62.67   Median :3.445   Median :182.6   Median : 7.874  
##  Mean   : 60.66   Mean   :3.359   Mean   :171.6   Mean   : 7.999  
##  3rd Qu.: 74.94   3rd Qu.:3.735   3rd Qu.:212.7   3rd Qu.: 8.622  
##  Max.   :100.14   Max.   :4.468   Max.   :227.2   Max.   :14.313  
## 
str(walmart) 
## spc_tbl_ [6,435 × 8] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Store       : Factor w/ 45 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Date        : Date[1:6435], format: "2010-02-05" "2010-02-12" ...
##  $ Weekly_Sales: num [1:6435] 1643691 1641957 1611968 1409728 1554807 ...
##  $ Holiday_Flag: Factor w/ 2 levels "0","1": 1 2 1 1 1 1 1 1 1 1 ...
##  $ Temperature : num [1:6435] 42.3 38.5 39.9 46.6 46.5 ...
##  $ Fuel_Price  : num [1:6435] 2.57 2.55 2.51 2.56 2.62 ...
##  $ CPI         : num [1:6435] 211 211 211 211 211 ...
##  $ Unemployment: num [1:6435] 8.11 8.11 8.11 8.11 8.11 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   Store = col_double(),
##   ..   Date = col_character(),
##   ..   Weekly_Sales = col_double(),
##   ..   Holiday_Flag = col_double(),
##   ..   Temperature = col_double(),
##   ..   Fuel_Price = col_double(),
##   ..   CPI = col_double(),
##   ..   Unemployment = col_double()
##   .. )
##  - attr(*, "problems")=<externalptr>

El dataset contiene 6.435 observaciones correspondientes a ventas semanales en cada una de las 45 tiendas Walmart en EEUU entre 2010 y 2012 además de 8 variables (de las cuales son 2 cualitativas y 6 cuantitativas).

A continuación, la descripción de cada una de las variables:

  • Store: Número de tienda que van del 1 al 45

  • Date: Día de la semana donde se recoge el dato de ventas. Tiene formato dd-mm-yyyy. Comienza 05-02-2010 hasta 26-10-2012; un total de 143 semanas (2 años y 9 meses)

  • Weekly_Sales: Ventas totales durante la semana

  • Holiday_Flag: Si la semana tiene un feriado especial o no. 1 si la semana tiene un feriado 0 si la semana es completamente laboral

  • Temperature: Temperatura media de la semana donde se producen las ventas

  • Fuel_Price: Precio del combustible en la región donde se encuentra situada la tienda

  • CPI: Índice de precios del cliente. Es un indicador utilizado para medir los cambios en el nivel de precios de bienes y servicios a lo largo del tiempo. El CPI es calculado por la Oficina de Estadísticas Laborales (Bureau of Labor Statistics) y es ampliamente utilizado como una medida de la inflación en la economía. Se basa en una “cesta de bienes” que representa los productos y servicios típicos que consume un hogar promedio. Es lo mismo que el IPC (Índice de Precio al Consumo).

  • Unemployment: Tasa de desempleo en la región donde se encuentra situada la tienda

15.2 Análisis y gráficos univariante

Comenzamos con unas tablas y gráficos simples por variable para ver como se comportan los datos y analizar los estadísticos principales por tienda.

# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
  format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}

# Calcular estadísticas de Walmart para la variable Temperature
walmart_stats_temperature <- walmart %>%
  group_by(Store) %>%
  summarise(min = formatCurrency(round(min(Temperature),2)),
            max = formatCurrency(round(max(Temperature),2)),
            median = formatCurrency(round(median(Temperature),2)),
            mean = formatCurrency(round(mean(Temperature),2)),
            sd = formatCurrency(round(sd(Temperature),2)))
            
# Calcular estadísticas de Walmart para la variable Fuel_Price
walmart_stats_fuel <- walmart %>%
  group_by(Store) %>%
  summarise(min = formatCurrency(round(min(Fuel_Price),2)),
            max = formatCurrency(round(max(Fuel_Price),2)),
            median = formatCurrency(round(median(Fuel_Price),2)),
            mean = formatCurrency(round(mean(Fuel_Price),2)),
            sd = formatCurrency(round(sd(Fuel_Price),2)))

# Calcular estadísticas de Walmart para la variable CPI
walmart_stats_cpi <- walmart %>%
  group_by(Store) %>%
  summarise(min = formatCurrency(round(min(CPI),2)),
            max = formatCurrency(round(max(CPI),2)),
            median = formatCurrency(round(median(CPI),2)),
            mean = formatCurrency(round(mean(CPI),2)),
            sd = formatCurrency(round(sd(CPI),2)))

# Calcular estadísticas de Walmart para la variable Unemployment
walmart_stats_unemployment <- walmart %>%
  group_by(Store) %>%
  summarise(min = formatCurrency(round(min(Unemployment),2)),
            max = formatCurrency(round(max(Unemployment),2)),
            median = formatCurrency(round(median(Unemployment),2)),
            mean = formatCurrency(round(mean(Unemployment),2)),
            sd = formatCurrency(round(sd(Unemployment),2)))

# Mostrar la tabla con los resultados para la variable Temperature
datatable(walmart_stats_temperature, options = list(scrollY = "300px", paging = FALSE), caption = "TEMPERATURE")
# Mostrar la tabla con los resultados para la variable Fuel_Price
datatable(walmart_stats_fuel, options = list(scrollY = "300px", paging = FALSE), caption = "FUEL PRICE")
# Mostrar la tabla con los resultados para la variable CPI
datatable(walmart_stats_cpi, options = list(scrollY = "300px", paging = FALSE), caption = "CPI")
# Mostrar la tabla con los resultados para la variable Unemployment
datatable(walmart_stats_unemployment, options = list(scrollY = "300px", paging = FALSE), caption = "UNEMPLOYMENT")
plot1 <- ggplot(walmart, aes(x = "", y = Temperature)) +
  geom_boxplot(fill = "#FFB21C") +
  labs(x = "ºF", y = "Temperature") +
  theme_minimal()

plot2 <- ggplot(walmart, aes(x = "", y = Fuel_Price)) +
  geom_boxplot(fill = "#FFB21C") +
  labs(x = "$ Gallon", y = "Fuel Price") +
  theme_minimal()

plot3 <- ggplot(walmart, aes(x = "", y = CPI)) +
  geom_boxplot(fill = "#FFB21C") +
  labs(x = "Index", y = "CPI") +
  theme_minimal()

plot4 <- ggplot(walmart, aes(x = "", y = Unemployment)) +
  geom_boxplot(fill = "#FFB21C") +
  labs(x = "Rate", y = "Unemployment") +
  theme_minimal()

grid <- grid.arrange(plot1, plot2, plot3, plot4, nrow = 2, ncol = 2)

Nos llama la atención la amplitud de las variables CPI y Unemployment y pasamos a graficarlas por tienda.

plot5 <- ggplot(walmart, aes(x = as.factor(Store), y = CPI, fill = as.factor(Store))) +
  geom_boxplot(fill = "#FFB21C") +
  labs(x = "Store", y = "CPI", title = "CPI per Store's Region") +
  theme_minimal() +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_flip()

plot6 <- ggplot(walmart, aes(x = as.factor(Store), y = Unemployment, fill = as.factor(Store))) +
  geom_boxplot(fill = "#FFB21C") +
  labs(x = "Store", y = "Unemployment", title = "Unemployment Rate per Store's Region") +
  theme_minimal() +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_flip()

plot5

plot6

Parece que dependiendo de la región donde se encuentre la tienda hay una mayor dispersión del gasto promedio en una cesta de la compra; y la tasa de desempleo también contribuye de manera relevante acentuándo dicha dispersión. El inconveniente es que a través del dataset desconocemos las regiones de EEUU donde se sitúan las 45 tiendas Walmart por lo que no podemos interpretar si tiene o no demasiado sentido, vamos a suponer que sí ya que no es lo mismo una tienda situada en un estado con mayor población que otro, y no es lo mismo una tienda en los estados rurales que en las grandes urbes de ambas costas.

Ahora nos disponemoa a analizar las ventas semanales usando la variable Weekly_Sales.

# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
  format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}

# Calcular estadísticas de Walmart
walmart_stats <- walmart %>%
  group_by(Store) %>%
  summarise(min = formatCurrency(round(min(Weekly_Sales),2)),
            max = formatCurrency(round(max(Weekly_Sales),2)),
            median = formatCurrency(round(median(Weekly_Sales),2)),
            mean = formatCurrency(round(mean(Weekly_Sales),2)),
            sd = formatCurrency(round(sd(Weekly_Sales),2)),
            sum = formatCurrency(round(sum(Weekly_Sales),2)))

# Mostrar la tabla con los resultados
datatable(walmart_stats, options = list(scrollY = "300px", paging = FALSE), caption = "SALES")
# Suma de las ventas semanales por tienda
ventas_por_tienda <- aggregate(Weekly_Sales ~ Store, data = walmart, FUN = sum)

# Orden de las tiendas de menor a mayor ventas
ventas_por_tienda <- ventas_por_tienda[order(ventas_por_tienda$Weekly_Sales, decreasing = FALSE), ]
ventas_por_tienda$Store <- factor(ventas_por_tienda$Store, levels = ventas_por_tienda$Store)

plot7 <- ggplot(ventas_por_tienda, aes(x = Weekly_Sales, y = Store)) +
  geom_bar(stat = "identity", fill = "#FFB21C") +
  labs(x = "Total Sales (Millions)", y = "Store", title = "Total Sales per Store Aggregated by Sum") +
  scale_y_discrete(labels = ventas_por_tienda$Store) +
  scale_x_continuous(labels = scales::unit_format(unit = "M", scale = 1e-6)) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

plot7

plot8 <- ggplot(walmart, aes(x = as.factor(Store), y = Weekly_Sales, fill = as.factor(Store))) +
  geom_boxplot(fill = "#FFB21C") +
  labs(x = "Store", y = "Weekly Sales", title = "Weekly Sales per Store") +
  theme_minimal() +
  theme(legend.position = "none") +
  coord_flip() +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6))

plot8

Este boxplot tampoco es muy exclarecedor, en muchas tiendas exite gran dispersión de los datos dentro del rango intercuartílico y hay datos muy atípicos. Esto puede sugerir que la distribución de datos está sesgada. Insistimos en que sin conocer las regiones donde se encuentran las tiendas es muy complicado llegar a más conclusiones o a un mejor entendimiento.

15.3 Análisis y gráficos multivariante

Una vez analizadas todas las variables de forma individual, podemos buscar algunas relaciones entre ellas mediante el siguiente análisis multivariante donde se pueden observar que existen correlaciones muy débiles, sólo pudiendo destacarse la correlación negativa del -0.3 que existe entre el CPI y la tasa de Unemployment; correlación que tiene bastante lógica dado que a mayor tasa de desempleo el gasto promedio en una cesta de la compra se verá redicido. Otro correlación muy débil aunque reseñable es la de 0.18 que existe entre la temperatura y el CPI, quizas a mayor temperatura la población salga más a la calle a consumir… pero es dificilmente demostrable.

variables <- c("Weekly_Sales", "Temperature", "Fuel_Price", "CPI", "Unemployment")
data <- walmart[, variables]
correlation_matrix <- cor(data)
plot9 <- ggcorrplot(correlation_matrix, type = "lower", lab = TRUE, lab_size = 3)

plot9

15.4 Análisis y gráficos de series temporales

Antes de graficar las series temporales de nuestro dataset vamos a sacar por tienda las 3 mayores ventas y el restos de las variables asociadas a las mismas para tratar de identificar si hay algún patrón como que todas se producen en las mismas fechas, que todas tienen unos índeces CPI concretos o que están sujetas a unas tasas de desempleo específicas…

# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
  format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}

# Calcular las 5 mayores cifras de ventas para cada tienda
walmart_top5sales <- walmart %>%
  group_by(Store) %>%
  top_n(3, Weekly_Sales) %>%
  summarise(Max_Weekly_Sales = formatCurrency(Weekly_Sales),
            Date = Date,
            Fuel_Price = formatCurrency(round(Fuel_Price,2)),
            CPI = formatCurrency(round(CPI,2)),
            Unemployment = formatCurrency(round(Unemployment,2)))

# Mostrar la tabla con los resultados
datatable(walmart_top5sales, options = list(scrollY = "300px", paging = FALSE), caption = "TOP 3 WEEKLY SALES PER STORE")

Ahora sacamos las 10 mayores ventas del dataset.

# Función para formatear números como moneda en formato americano
formatCurrency <- function(x) {
  format(x, big.mark = ".", decimal.mark = ",", nsmall = 2, trim = TRUE, scientific = FALSE)
}

# Obtener los 10 primeros días con más ventas y aplicar el formato de moneda
top_10_days <- walmart %>%
  group_by(Date) %>%
  summarise(total_sales = sum(Weekly_Sales)) %>%
  top_n(10, total_sales) %>%
  arrange(desc(total_sales)) %>%
  mutate(total_sales = formatCurrency(total_sales))

datatable(top_10_days, options = list(scrollY = "300px", paging = FALSE), caption = "TOP 10 DAYS WITH HIGHEST SALES")

Se observa claramente como en ambas tablas las semanas de Navidad y Acción de Gracias en noviembre tienen el record de ventas todos los años.


GRÁFICO SERIE TEMPORAL

Comenzamos realizando con unos gráficos de líneas temporales en las variables más relevantes que deseamos analizar.

walmart <- arrange(walmart, Date)

plot10 <- ggplot(walmart, aes(x = Date, y = Weekly_Sales)) +
  geom_line() +
  labs(x = "Date", y = "Weekly Sales", title = "Weekly Sales Time Series Plot") +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6)) + 
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

plot10

plot11 <- ggplot(walmart, aes(x = Date, y = Temperature)) +
  geom_line() +
  labs(x = "Date", y = "Temperature", title = "Temperature Time Series Plot") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_minimal()

plot12 <- ggplot(walmart, aes(x = Date, y = Fuel_Price)) +
  geom_line() +
  labs(x = "Date", y = "Fuel Price", title = "Fuel Price Time Series Plot") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_minimal()

plot13 <- ggplot(walmart, aes(x = Date, y = CPI)) +
  geom_line() +
  labs(x = "Date", y = "CPI", title = "CPI Time Series Plot") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_minimal()

plot14 <- ggplot(walmart, aes(x = Date, y = Unemployment)) +
  geom_line() +
  labs(x = "Date", y = "Unemployment", title = "Unemployment Time Series Plot") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme_minimal()

grid <- grid.arrange(plot11, plot12, plot13, plot14, nrow = 2, ncol = 2)

El gráfico de la temperatura muestra un comportamiento totalmente normal donde se pueden identificar a simple vista las estaciones y sus variaciones de temperatura entre invierno y verano.

El gráfico del precio del combustible muestra una relativa estabilidad en precios sufriendo una subida a partir del primer trimestre del 2010 teniendo su punto más álgido en el segundo trimestre del año; punto a partir del cual vuelve a bajar situándose cerca de los 2,5 dólares por galón en el tercer trimestre del año. No es así en el cuarto trimestre donde se origina una vertiginosa subida del precio por galón hasta el segundo trimestre del 2011, pasando de 2,5 a 4 dólares. Desde este punto hasta el final del año se produce otra severa caída de precios cerrando el año 2011 alrededor de los 3,25 dólares. Durante el año 2012 se van produciendo unas subidas y bajadas de precios con bastante apuntamiento entre los 3,25 y 4 dólares.

El gráfico del CPI muestra un comportamiento monótono creciente pero de manera apaisada que nos indica que el gasto medio en la lista de la compra va aumentando lentamente de manera progresiva. No obstante el gráfico muestra gran varianza en este índice debido precisamente a lo ya comentado con anteriorida; hay tiendas Walmart situadas en regiones de EEUU donde índice de población y el poder adquisitivo es menor y otras donde es mayor, diferencias entre estados rurales y grandes urbes.

El gráfico de la tasa Unemployment parece que se muestra estable a lo largo del 2010 y 2011 hasta comenzar 2012 donde generalmente se va produciendo una caída en las tasas de desempleo.

15.4.1 Estacionariedad

Antes de analizar la serie mediante test estadísticos, analicemos gráficamente más en detalle, para confirmar la estacionalidad detectada. Dadas las características de la serie realizamos una descomposición STL para corroborar que existe tendencia y estacionalidad.


La descomposición STL (Seasonal and Trend decomposition using Loess) es un método utilizado para descomponer una serie temporal en sus componentes de tendencia, estacionalidad y residuos. Utiliza un enfoque no paramétrico basado en el ajuste de suavizadores locales (Loess) para estimar la tendencia y la estacionalidad de la serie temporal. La tendencia representa el patrón general y de largo plazo en los datos, mientras que la estacionalidad captura los patrones recurrentes y periódicos. Los residuos son la parte de la serie temporal que no puede ser explicada por la tendencia y la estacionalidad.

# Convertir el dataset a un tsibble
walmart_tsibble <- as_tsibble(walmart, key = Store, index = Date)
walmart_tsibble_train <- walmart_tsibble %>% filter_index(. ~ "2012-06-29")
walmart_tsibble_test <- walmart_tsibble %>% filter_index("2012-07-06" ~ .)
walmart_tsibble_train %>%
  model(seats = feasts:::STL(Weekly_Sales)) %>%
  components() %>%
  autoplot() +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(labels = unit_format(unit = "M", scale = 1e-6)) +
  theme(panel.grid.major = element_line(color = "grey"),
        panel.grid.minor = element_line(color = "grey"),
        panel.border = element_rect(color = "grey", fill = NA))

Se confirman visualmente las estacionalidades supuestas.

walmart_tsibble_train %>%
  model(seats = feasts:::STL(Fuel_Price)) %>%
  components() %>%
  autoplot() +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.grid.major = element_line(color = "grey"),
        panel.grid.minor = element_line(color = "grey"),
        panel.border = element_rect(color = "grey", fill = NA))

walmart_tsibble_train %>%
  model(seats = feasts:::STL(CPI)) %>%
  components() %>%
  autoplot() +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.grid.major = element_line(color = "grey"),
        panel.grid.minor = element_line(color = "grey"),
        panel.border = element_rect(color = "grey", fill = NA))

walmart_tsibble_train %>%
  model(seats = feasts:::STL(Unemployment)) %>%
  components() %>%
  autoplot() +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.grid.major = element_line(color = "grey"),
        panel.grid.minor = element_line(color = "grey"),
        panel.border = element_rect(color = "grey", fill = NA))


Estacionariedad en MEDIA

walmart_tsibble %>%  features(Weekly_Sales, unitroot_kpss)
## # A tibble: 45 × 3
##    Store kpss_stat kpss_pvalue
##    <fct>     <dbl>       <dbl>
##  1 1        0.476       0.0471
##  2 2        0.0787      0.1   
##  3 3        0.608       0.0219
##  4 4        0.920       0.01  
##  5 5        0.672       0.0161
##  6 6        0.0365      0.1   
##  7 7        0.512       0.0390
##  8 8        0.187       0.1   
##  9 9        0.669       0.0163
## 10 10       0.131       0.1   
## # … with 35 more rows

Estacionariedad en VARIANZA

15.4.2 Determinación del modelo

15.4.2.1 Modelo A

15.4.2.2 Modelo B

15.4.3 Fase de estimación y contraste

15.4.4 Diagnosis de residuos

Test de Media = 0

Test de Autocorrelaciones

Test de Homocedasticidad

Test de Normalidad

15.4.5 Fase de predicción

Finalmente evaluamos la capacidad predictiva del modelo. Calculamos las predicciones y calculamos sus errores y los comparamos con los residuos.